Bayesian inversion of a diffusion model 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 modelling). When complex dynamical systems are considered, such as partial differential equations, this task may become challenging or ill-posed. In this work, a linear parabolic equation is considered as a model for protein transcription from MRNA. The objective is to estimate jointly the differential operator coefficients, namely the rates of diffusion and self-regulation, as well as a functional source. The recent Bayesian methodology for infinite dimensional inverse problems is applied, providing a unique posterior distribution on the parameter space continuous in the data. This posterior is then summarized using a Maximum a Posteriori estimator. Finally, the theoretical solution is illustrated using a state-of-the-art MCMC algorithm adapted to this non-Gaussian setting.


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 equation (with constant diffusion and damping without transport): where the spatial domain is an open set Ω ⊂ 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 z (hereafter called the solution of Eq. 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-defined solution [well-posedness in the sense of Hadamard Hadamard (1902)] but only a particular case will be dealt with here. 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 and the latter is the objective of this paper. The motivation comes from a challenging identification problem in Biology where the objective is to infer jointly self-regulation and diffusion rates with the source f given a limited number of noisy observations of the solution z. Note that given their physical interpretation, the parameters u = (λ, D, f ) must all be non-negative (in an obvious sense). This problem has already been solved using various approaches. In Becker et al. (2013), the authors use a system of ordinary differential equations instead of Eq. 1, and minimize a discrete version of a least-square penalty functional, while confidence intervals on parameters are given by bootstrapping. In a Bayesian setting, alternative methods have been based on Latent Force Models (Alvarez et al. 2013;Särkkä et al. 2018). These approaches assume that the source can be modelled with Gaussian Processes (Rasmussen and Williams 2006). In particular, if f is taken to be such a process and if the decay and the diffusion are constants, z is Gaussian as well (since z depends linearly in f ). The two constant parameters λ and D can then be optimized as hyper-parameters through standard likelihood maximization (Lopez-Lopera et al. 2019). Major drawbacks of this approach are both the difficulty to enforce positiveness of the reconstructed source and the absence of uncertainty quantification on λ, D.
In this work, a more general methodology is applied, based on the recent advent of Bayesian Inverse Problems (Stuart 2010) for infinite dimensional spaces. This has the advantage of dealing with the ill-posedness while fully integrating the quantification of uncertainties in a unified approach. Moreover, the possibility to include previous physical constraints in the prior will be particularly meaningful. The rest of this paper is organized as follows: Sect. 2 presents all the mathematical analysis underlying the Bayesian inversion, namely sufficient regularity of the forward operator (mapping the parameter (λ, D, f ) to the PDE solution z), existence and uniqueness of the posterior under a simple class of priors, and finally uniqueness of the associated posterior. Section 3 focuses on an adaptation of a geometric Markov Chain Monte-Carlo algorithm, well-defined in function spaces. Finally, Sect. 4 contains all implementation details, as well as a methodology to tune hyper-parameters of the model. All numerical results are based on a real-world dataset related to the developmental biology of the Drosophila Melanogaster.

Bayesian inversion
As previously announced, the 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 z (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 more constraints with Tikhonov-Philips functionals, to ensure uniqueness and continuity w.r.t. the observations (Isakov 2017;Schuster et al. 2012). Doing so, the regularized solution will be compatible with the dataset. 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 (Ghanem et al. 2017;Biegler et al. 2011;Sullivan 2015). One approach consists in treating these 2 objectives sequentially, first regularizing then quantifying the resulting uncertainty. However, the Bayesian methodology for inverse problems (Stuart 2010) and more recently Dashti and Stuart (2015) 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, which leads to robust algorithms w.r.t. discretization dimension. Indeed, finite approximations of probability measures may be absolutely continuous while their infinite counterparts are mutually singular. This becomes particularly troublesome in MCMC sampling for instance (Cotter et al. 2013).
In essence, instead of searching for one particular parameter solving the regularized problem, this approach considers conditional probability measures of parameters given observations. Namely, given a prior distribution (see Sect. 2.2) and few technical conditions on the forward operator (see Sect. 2.1), Bayes theorem applies and exhibits a unique posterior distribution (see Sect. 2.3), which is continuous in the data [w.r.t. Hellinger metric (Stuart 2010)]. Finally, one may summarize information from this posterior distribution with point estimators such as expected value or modes (Sect. 2.4). All these steps will now be presented.

Forward model analysis
The first step is to detail precisely the required regularity of the solution map from Eq. 1. Using common variational techniques from PDE theory [see Evans (1998) or Brezis (2011)], one can show that this equation has a unique weak solution (see Proposition 1 which proof is in the appendix) given u = (λ, D, f ) in a domain U that will be explicited later on. Moreover, this solution evolves smoothly when the parameter varies. Without any loss of generality and keeping in mind the biological application, the underlying physical domain will be Ω =]0, L[ with L ∈ R + , representing the anterior-posterior axis of a Drosophilia embryo.
Moreover, this map has the following properties: 1. (Energy estimate): it satisfies the following estimate ∀u ∈ U , where C > 0 is a constant independent of u, 2. (Local Lipschitz continuity): ∀u in the interior of U , ∀r > 0 such that B(u, r ) ⊂ U , ∃L(u, r ) > 0, The proof is standard using methods from PDE and optimal control theory (Evans 1998). The properties stated in proposition 1 will be important for subsequent analysis (Sect. 2.3): 1. The energy estimate is fundamental to establish most of the results concerning the posterior distribution. Indeed, it gives a precise upper bound, useful to show integrability of the yet-to-come likelihood, 2. Local Lipschitz continuity implies measurability of the solution map w.r.t. the Borel σ -algebra and is used in the characterization of posterior modes (Maximum a Posteriori estimators), 3. Second order Fréchet differentiability will be necessary for geometric methods in optimization (research of posterior modes) and sampling (Markov-Chain Monte-Carlo) since they rely on Hessian-type information.

Choice of a prior distribution
The second step is to choose a prior probability distribution on U , encoding all knowledge on the problem at hand, while being simple enough to keep the analysis tractable. The prior will be constructed as a product measure, specifying marginal measures μ λ 0 , μ D 0 , μ f 0 on all parameters: For simplicity, both measures for λ and D will be taken uniform on respective intervals . Now, since f must be non-negative, 1 is re-parametrized with the following source term: where In this work, μ f 0 will be taken as the Gaussian measure associated with a continuous Gaussian process [see Bogachev (1998)

for a presentation of infinite dimensional
Gaussian measures] such that f is almost-surely continuous. Remark that the exponential map in Eq. 3 could be replaced with any sufficiently differentiable function from R to R + (thus keeping the second order Fréchet differentiability of the solution map). Besides, alternative distributions are also possible for f like Besov priors [from Dashti et al. (2013)] or more general convex measures [from Hosseini and Nigam (2017)]. However, this choice is also motivated by practical reasons, since one can build a Gaussian measure μ re f dominating μ 0 : The parameters of μ re f are tuned by choosing (minimizing Kullback-Leibler divergence of μ re f relative to μ 0 ). This dominant measure will be critical for posterior modes analysis (Sect. 2.4) and MCMC sampling (Sect. 3.2).

Posterior distribution
It is now time to show that the particular setting provided so far (forward model and prior distribution) leads to a well defined posterior measure. This is the purpose of Proposition 2 which is a direct application of the theory initially developed in Stuart (2010) [see Dashti and Stuart (2015) for an updated presentation]. In this purpose, consider a dataset y = (y i ) i∈ [1,n] ∈ R n consisting of observations from the solution z at different times and locations (t i , x i ) i∈ [1,n] ∈ ([0, T ] × [0, L]) n , produced under the following additive noise model (in vector notations): 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 (z[u](x i , t i )) i∈ [1,n] (composition of solution map z with Dirac type measure). Note that it is well-defined since the function z is continuous on the domain for every u ∈ U . The following proposition, which is again proved in the appendix, establishes the existence, uniqueness and continuity in y of the posterior probability measure μ y (the solution of the inverse problem), expressing how observations y updated prior beliefs μ 0 on the parameter u.
Proposition 2 Let G be the observation operator defined in Eq. 4, y ∈ R n a dataset and μ 0 the probability measure defined in Eq. 2, then there exists a unique posterior measure μ y for u|y. It is characterized by the following Radon-Nikodym density w.r.t. μ 0 : with the negative log-likelihood and the marginal Furthermore, μ y is continuous in y w.r.t. Hellinger distance.
In fact, Proposition 2 gives two distinct results: a) the existence and uniqueness of a posterior (as long as μ 0 is Radon and μ 0 (U ) = 1, which is the case here), b) wellposedness of the Bayesian inverse problem. In particular, the use of a Gaussian prior on f gives sufficient integrability, even under the re-parametrization from Eq. 3 (using Fernique's theorem). If one chooses a different map between f * and f in Eq. 3, this condition may be considerably relaxed (using something slower than the exponential) and prior measures with lower integrability conditions can be considered.

Maximum a posteriori estimator
In the previous section, the well-posedness of the Bayesian inverse problem has been proved. However, the posterior distribution is only known up to a multiplicative constant, through its density w.r.t. μ 0 . In the application, μ y will need to be summarized, which is usually done by the selection of a particular estimator in U , such the posterior mean or a mode. One consequence of Proposition 2 is that the posterior mean is automatically continuous in the data y [since well-posedness is w.r.t. Hellinger distance, see Dashti and Stuart (2015)]. However, optimality properties (in a decision theoretic context) are not yet well-understood in infinite dimension to the best of our knowledge. This is why posterior modes (or Maximum a Posteriori) are considered instead. Furthermore, they provide a clear link with the classical Tikhonov-Philips regularization [see Dashti et al. (2013); Helin and Burger (2015); Agapiou et al. (2018)] and a useful variational characterization (in case of Gaussian or Besov priors) which will be the cornerstone of the numerical application, see Proposition 3 (which proof is given in the appendix).
Proposition 3 Let μ 0 be the prior probability measure defined in Eq. 2 and μ re f the Gaussian reference measure from Eq. 2.2, then the modes of μ y are exactly the minimizers of the following (generalized) Onsager-Machlup functional: where . μ f 0 is the Cameron-Martin norm associated to μ f 0 .
A minimizer of the previous generalized Onsager-Machlup functional will be noted u M AP = (λ M AP , D M AP , f M AP ) and need not be unique. The precise application of this proposition to the biological setting is done in Sect. 4.1, once the prior distribution is fully specified.

Approximation
The final step in the theoretical analysis of the Bayesian inverse problem is to study its approximation properties, since it will be solved numerically. There are two important things to check: (a) properties of the approximated posteriors (since it is what is available), (b) the consistency of these approximated posteriors. This will be done by projection of the parameter u onto a finite dimensional subspace of U , constructed with a stochastic basis under μ f 0 [see Okazaki (1986) for a detailed introduction on stochastic bases including Banach spaces]. The second source of approximation is the use of a numerical solver for the PDE solution, but this will be neglected here (but can be conducted in a subsequent work). The next proposition will establish that the posterior distribution is well approximated, giving an estimate of the error when the Hilbert basis is chosen specifically (spectral basis of the covariance operator, namely Karhunen-Loève decomposition) and its proof is in the appendix.

Metropolis-Hastings algorithm
As it was previously announced, the main motivation for the Bayesian methodology here is quantification of uncertainty, which will be done by sampling the posterior measure. Among the vast catalogue of methods for probability distributions simulation (Sequential Monte-Carlo, Approximate Bayesian Computations, Transport Maps, etc...), Markov chain Monte-Carlo is very popular [MCMC, see Brooks et al. (2011)] and well defined on function spaces Tierney (1998) even though ergodicity analysis of such algorithms is still in its infancy (Hairer et al. 2005(Hairer et al. , 2007(Hairer et al. , 2014Rudolf and Sprungk 2018). After a short presentation of the Metropolis-Hastings algorithm (Sect. 3.1), this section will focus on a state-of-the-art Markov kernel designed to sample from Gaussian measures (Sect. 3.2) and adapt it to the current non-Gaussian prior using the Gaussian dominating measure μ re f .

Metropolis-Hastings on function spaces
The Metropolis-Hastings algorithm (MH) is a very general (Tierney 1998) method to design Markov chains to sample from a given probability measure. It is based on a two-step process on each iteration: 1. 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 ), 2. Accept the new state v with probability α (u, v) or remain at u.
This algorithm provides a sample distributed under a predefined probability measure μ, if one selects α and Q in a specific way [see Tierney (1998) for a discussion in general state spaces]. For instance, let ν(du, dv) = μ(du)Q(u, dv) and ν T (du, dv) = ν(dv, du), the Metropolis-Hastings algorithm typically considers the following acceptance probability: which, in particular, requires the absolute continuity of ν T w.r.t. ν (detailed balance condition 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 Dashti and Stuart (2015); Girolami and Calderhead (2011); Beskos et al. (2017); Cotter et al. (2013); Hairer et al. (2014)] is to select Q revertible w.r.t. μ 0 . Indeed, in this case (with ν 0 (du, dv) = μ 0 (du)Q(u, dv)): In theory, the MH algorithm can 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: -to adjust to locally mimic the target distribution μ y , -include a step size to tune acceptance probability to reasonable values.
These two properties may be used to trade-off self-correlation, acceptance rates and convergence speed to high interest areas of the parameter space. Next section presents a proposal Q with both properties to sample from a Gaussian prior distribution.

Geometric MCMC under Gaussian reference
The specific Markov proposal kernel Q, tailored to sample distributions having a density w.r.t. a Gaussian measure μ re f will now be presented. Most of the recent work on infinite dimensional MCMC methods is based on the following Langevin stochastic differential equation: 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 Beskos et al. (2017), a semi-implicit discretization of Eq. 7 leads to a Markov chain with the following kernel: where h > 0 is a step-size parameter, ρ = 1−h 1+h and: This dynamic explores the parameter space with a balance between Newton-type 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 [Law (2014) and later generalized in Beskos et al. (2017) and Cui et al. (2016)]. 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 for instance) become ineffective in this case. However, the infinite dimensional manifold Modified Adjusted Langevin Algorithm [∞-mMALA from Beskos et al. (2017)] considers a specific preconditioner: where H Φ (u) is the Gauss-Newton Hessian operator of Φ, which locally adapts to the posterior. This kernel does not preserve the distribution μ y but is shown to be absolutely continuous w.r.t. the reference measure μ re f , almost-surely in u (under technical assumptions regarding K (u) linked with Feldman-Hajek theorem) and the Radon-Nikodym density is: as it is done in Beskos et al. (2017), it finally comes: Finally, the acceptance probability associated to the Markov kernel from Eq. 8 is 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 Beskos et al. (2017) or likelihood-informed in Cui et al. (2016)] to reduce the computational burden. A second alternative is to choose a constant preconditioner, located at a (precomputed) posterior mode for instance [similar to HMALA in Cui et al. (2016) and gpCN in Rudolf and Sprungk (2018)], which is done in this work.

Numerical application
This section will now introduce the practical implementation of the previous methodology on the problem of reverse-engineering for post-transcriptional gap-gene in Drosophila Melanogaster. First, the prior distribution is detailed, as well as a random series representation for f , leading to an approximation of f * . Moreover, the associated generalized Onsager-Machlup functional is also provided. Then, quantitative results are given on the dataset taken from Becker et al. (2013), consisting in protein concentration measurements irregularly spread in space and time.

Choice of a continuous Gaussian process
In the previous analysis, μ f 0 has been defined as the probability measure over C ([0, L] × [0, T ]; R) associated with a continuous Gaussian process. In practice, it will be chosen centred for simplicity, thus completely specified by a covariance kernel over the space-time domain [0, L] × [0, T ]. The literature on such processes is vast [see Rasmussen and Williams (2006) for instance], and here a tensor product of two Brownian bridges (in time and space) will be used. The associated covariance kernel is given by: where the scaling constant is such that the process has maximum variance equal to σ 2 . In particular, this process has a well-known Karhunen-Loeve decomposition (see Bay and Croix) using Schauder-type hat functions. In turn, the process will thus be approximated as follows: where (ξ i 1 ,i 2 ) 1≤i 1 ,i 2 ≤N are i.i.d. N (0, 1) random variables and ϕ i are hat functions on dyadic intervals and the precise weights are given in Bay and Croix. We consider N = 20 (400 basis functions) thus the approximated parameterũ = (λ, D,f ) is of dimension 402.

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 Alnaes et al. (2015) and Langtangen and Logg (2017)] and finite differences in time. We use 100 finite elements and 30 time steps on a 2018 standard laptop computer. 1 We set L = 100 and T = 100 is the final time (data-points have been rescaled). All quantities related to negative log-likelihood derivatives (Gradient and Gauss-Newton Hessian matrix) are numerically computed using discrete adjoint methods [see Hinze et al. (2009) or Heinkenschloss (2008] to keep scalability in N . The initial point in the chain is chosen at the MAP location, obtained by minimization of the functional in Eq. 3 (Prior based initialization results in long burnin phase). Practical optimization is done using L-BFGS-B algorithm from the Scipy library (Byrd et al. 1995).

Estimation of noise variance
So far, the prior measure as well as the observation model include hyper-parameters playing fundamental roles and needing to be tuned. The complete list is λ M , D m , D M , σ 2 η and to the best of our knowledge, there are no general methods to estimate them efficiently in this context. Indeed, no closed formulae exists for the likelihood (probability density of G (u)) and cross-validation seems computationally out of reach. However, we provide here an empirical approach to tune the noise variance level σ 2 η . The other parameters are fixed to arbitrary values. Our approach consists in: 1. Fixσ 2 η = 1. 2. Find u M AP minimizing I using the current noise level σ 2 η . 3. Update the current variance estimationσ 2 η = 1 n−1 y − G (u M AP ) 2 . 4. Go back to 2. The algorithm is stopped once the parameter σ 2 η reaches a stable value, which in this application isσ 2 η = 20, 50. The estimator u M AP is taken as the last computed minimizer of the (discretized) Onsager-Machlup functional I .

Results
We now turn to our main objective, the inversion and uncertainty quantification of gap-gene protein concentration from Becker et al. (2013). The dataset consists of 508 different measures which are non-uniformly spread in time and space (precise repartition can be seen in Fig. 1). 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 21,000 total iterations and the resulting traceplot is given in Fig. 4 for negative log-likelihood, decay, diffusion and first three components of f . The first thousand iterations are used as burnin and according to the autocorrelation function (Fig. 3), we choose to keep one iteration out of two hundreds as posterior sample (thinning). From this, we compute both posterior mean and MAP estimates, the precise values of decay, diffusion, negative log-likelihood and Onsager-Machlup functional being given in Table 1.  Additionally to the estimated values, one can also look at the marginal distribution on Fig. 2. These values suggest that there is lack of identificability between Λ and D (related to the non-injectivity of the forward operator). However, the strength of the Bayesian approach is to quantify this phenomenon, giving here potential values compatible with the data. Concerning the MAP estimator (Fig. 1), we recover both events described in Becker et al. (2013), 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 Fig. 5). However, the point-wise variance on the solution y remains important before the first observations, which translates into an expected level of uncertainty. This indicates that given the data (and more precisely its time/space repartition), multiple scenarios are valid and our model would require more data to progress in these areas.

Conclusions
In this work, we applied the Bayesian inverse problem methodology from Stuart (2010) 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 MAP estimates, a consistent approximation and apply a state-of-the-art MCMC methodology. Because the forward MAP is non-linear, the uniqueness 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 precious 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 intend to address it in a future work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Proofs
Proof of proposition 1 Using standard notations from PDE theory, let equipped with the norm: 1. (Existence, uniqueness and energy estimate) with λ M > 0 and 0 < D m ≤ D M < ∞, then the differential operator in Eq. 1 is uniformly parabolic. By standard techniques, e.g. Galerkin approximation, one can show the existence and uniqueness of a weak solution w ∈ W to Eq. 1 [Theorems 3 and 4, section 7.1, Evans (1998)]. Moreover, the regularity of the source term implies an improved regularity for the weak solution, namely z ∈ L 2 ([0, T ], H 2 (Ω)) with z ∈ L 2 ([0, T ], H ) [Theorem 5, section 7.1, Evans (1998)]. The last result is an energy estimate: Since D is bounded below by D m > 0, the constant C only on the coercivity property of the differential operator, thus D m and does not depend on D. This in turn implies the announced energy estimate, by the two continuous embeddings: [Theorem 4, section 5.9, Evans (1998)]: 2. (Local Lipschitz continuity in u, weaker norm) Let u in the interior of U , r > 0 such that B(u, r ) ⊂ U and (u 1 , u 2 ) ∈ B(u, r )×B(u, r ). There exists two unique solutions with respect to u 1 and u 2 such that ∀i ∈ {1, 2}, ∀v ∈ L 2 ([0, T ], V ) and for almost-every t in [0, T ]: which, using subtraction and rearranging the terms, leads to: Now, let v = z 1 − z 2 in Eq. 14, drop the term λ 1 z 1 − z 2 2 H in the left-hand side, use Cauchy-Schwarz and Poincaré's inequalities in the right-hand side, then: where c ≥ 0 is Poincaré's constant. Integrating both sides between 0 and T , using the identity d Now, it remains to use the identity ab ≤ a 2 2 2 + 2 b 2 2 , that is with C(u, r ) ≥ 0. Let's rewrite equality 14 this way: Applying Cauchy-Schwarz and Poincaré's inequalities in the right-hand side of Eq. 17 gives: Using integration and Cauchy-Schwarz inequality, one obtains which concludes the proof on the local Lipschitz continuity. Finally, the continuous embedding of C ([0, L] × [0, T ], R) into the more regular space, gives the result. 4. (Fréchet differentiability) Now, we deal with the Fréchet differentiability. Let u = (λ, D, f ) in the interior of U , h u = (h λ , h D , h f ) ∈ U such that u + h u ∈ U and h z ∈ W , then ∀v ∈ L 2 ([0, T ], V ), we will use the following notation: then with (using triangle, Poincaré's and Cauchy-Schwarz inequalities): where C ≥ 0 is a (new) constant independent from u and: Using again Cauchy-Schwarz and Poincaré's inequalities, we have: with C ≥ 0 another constant independent of (h z , h u ) thus F (z, u) is linear and bounded, which means that F is Fréchet-differentiable on W ×U . Consider F z the This now gives ∀y ∈ B(0, r ): Now, one can use the linearity of z in f * , thus It remains to see that and since f − P m f E → 0, the result is clear as Fernique's theorem in E gives the required integrability.