The Inverse Problem of Reconstructing Reaction-Diffusion Systems

This paper considers the inverse problem of recovering state-dependent source terms in a reaction-diffusion system from overposed data consisting of the values of the state variables either at a fixed finite time (census-type data) or a time trace of their values at a fixed point on the boundary of the spatial domain. We show both uniqueness results and the convergence of an iteration scheme designed to recover these sources. This leads to a reconstructive method and we shall demonstrate its effectiveness by several illustrative examples.


Introduction
Reaction diffusion equations have a rich history in the building of mathematical models for physical processes. They are descendants of nonlinear ordinary differential equations in time with an added spatial component making for a partial differential equation of parabolic type. These early models dating from the first decades of the twentieth century include that of Fisher in considering the Verhulst logistic equation together with a spatial diffusion or migration, u t − ku xx = f (u) = bu(1 − cu) to take into account migration of species and that of Kolmogorov, Petrovskii and Piskunov in similar models which are now collectively referred to as the Fisher-KPP theory of population modeling, see, [21]. There is also work in combustion theory due to Zeldovich and Frank-Kamenetskii that utilize higher order polynomials in the state variable u and where the diffusion term acts as a balance to the chemical reactions, [8].
The use of systems of reaction diffusion models followed quickly; adding a spatial component to traditional population dynamic models such as predator-prey and competitive species as well as the interaction of multiple species or chemicals. By the early 1950's it was recognized by Alan Turing that solutions to such equations can, under the correct balance of terms, be used to simulate natural pattern formations such as stripes and spots that may arise naturally out of a homogeneous, uniform state [30]. This theory, which can be called a reaction-diffusion theory of morphogenesis, has been a major recurrent theme across many application areas.
These models use the underlying physics to infer assumptions about the specific form of the reaction term f (u). The few constants appearing, if not exactly known, are easily determined in a straightforward way by a least squares fit to data measurements. We envision a more complex situation where the function f (u) (or multiple such functions as we will be considering systems of equations) cannot be assumed to have a specific known form, or to be analytic so that knowing it over a limited range gives a global extension, and therefore must be treated as an undetermined coefficient problem for a nonlinear partial differential equation.
Let Ω be a bounded, simply-connected region in R d with smooth boundary ∂ Ω, and let L := −∇ · (a(x)∇·) + q(x). be a uniformly elliptic operator of second order with L ∞ coef- for some fixed time T and subject to the prescribed initial and boundary conditions ∂ u ∂ ν (x,t) + γ 1 (x)u(x,t) = θ 1 (x,t) where ν is the outer unit normal and γ i (x), θ i (x,t) are in C β (∂ Ω), where C β is the Schauder space with Hölder exponent β , 0 < β < 1.
In equation (1) we assume r 1 (x,t, u, v) and r 2 (x,t, u, v) are known and the interaction variable w = w(u, v) is also known but either the pair { f 1 , f 2 } or the pair {φ 1 , φ 2 } is unknown and the inverse problems posed are to determine these quantities.
Thus, there are two distinct inverse problems. The first is when we assume the interaction coupling φ i (w) between u and v is known, but both f 1 (u) and f 2 (v) have to be determined.
The second is when we assume the growth rate couplings f i for both u and v are known, but the interaction terms φ 1 (w) and φ 2 (w) have to be determined.
In order to perform these recoveries we must prescribe additional data and we shall again consider two possibilities: the values of u(x, T ) and v(x, T ) taken at a later, fixed time T u(x, T ) = g u (x), v(x, T ) = g v (x), x ∈ Ω; or the time traces u(x 0 ,t) and v(x 0 ,t) measured at a fixed point x 0 ∈ ∂ Ω and for all t ∈ (0, T ).
If the boundary conditions at the measurement point x 0 are of Dirichlet type (γ = ∞) then instead we measure the flux a(x 0 )∂ ν at x 0 in (4). Note that final time data corresponds to census data taken a fixed time T for the species involved. The time trace data involves monitoring the population (or of chemical concentrations) at a fixed spatial point as a function of time. Both of these data measurements are quite standard in applications. In (3), the observation domain can be restricted to a subdomain ω of Ω; in view of the fact that the functions to be discovered are univariate, even an appropriately chosen curve in Ω could possibly suffice.
The interaction coupling w of u and v, which we assume known, can take on several forms. The near universal choice in ecological modeling is to take w = uv. This is also common in other applications but other more complex possibilities are in use.
For example, the Gray-Scott model of reaction diffusion, [23], takes w = u 2 v and is a coupling term that often leads to pattern formation. This coupling occurs, for example, in molecular realizations where there is an activator (u) and an inhibitor (v). The antagonistic effect here occurs from the relative depletion of v that is consumed during the production of u. The so-called Brusselator equation (the Walgraef and Aifantis equation) which also leads to the generation of sustained oscillations, instead takes w = bu − cu 2 v and occurs, amongst many other situations, in the dislocation dynamics in materials subjected to cyclic loading, [31].
Other possibilities include w = √ u 2 + v 2 or the nonlocal situation w = Ω (u 2 + v 2 ) dx as well as combinations of all of the above.
The astute reader will already have noted that there are other very similar combinations that might lead to physically-motivated inverse problems. Indeed, this is the case and the analysis from Section 3.1 also applies to other combinations of unknowns than those mentioned above.
In the case of a single equation using time trace data, uniqueness results and the convergence of reconstruction algorithms were shown in [4,27,25] for the recovery of the unknown term f (u). In the more recent paper [16], the authors used final time data as in (4) and showed uniqueness and contractibility results for a fixed point iteration scheme which allowed effective recovery of f (u). We also point to [15] where the reaction term is of the form q(x) f (u) with known f (u) and an unknown space dependent term q(x) that controls the intensity of the reaction.
The reconstruction methods that we consider here are based on a projection of the PDE on the observation manifold, which naturally leads to a fixed point iteration for the unknown terms. As compared to Newton's method, which is a more general approach, this projection principle requires an appropriate form of the given data; in particular, the observation manifold cannot be orthogonal to the dependence direction of the unknown function, i.e., if we aim at reconstructing an x dependent coefficient then time trace observations cannot be used in such a projected approach (and will hardly give good reconstruction results with any other method either). Since here the unknown function does not depend on x or t this problem does not appear, though. Moreover, a convergence proof of regularized Newton type methods can only be carried out under certain conditions on the forward operator, such as the tangential cone condition (see, e.g., [14] and the references therein), which is most probably not satisfied here. Also note that the fixed point approach used here entails a uniqueness result for the inverse problem.
We remark that we have followed standard mathematical practice when writing differential equations by scaling the equation to set parameters not under consideration to unity. Thus in our reconstruction we have taken the volume of Ω to be one and the elliptic operator L to be the Laplacian thereby making the implicit assumption that the diffusion coefficient a(x) in the leading term div(a∇u) is unity. In an actual physical situation this will not be the case; for example the value of a for a molecule diffusing in a gas will be several orders of magnitude smaller.
The outline of the paper is as follows. Section 3 provides an extension of the convergence analysis in [16] for a fixed point iteration as well as its consequences for uniqueness in two directions. Firstly, we "vectorize" this analysis to be applicable for a class of reactiondiffusion equations, characterized by certain conditions. Secondly, we allow for higher space dimensions by carrying out the estimates in Schauder spaces rather than (Hilbert) Sobolev spaces, thus largely avoiding dependence on the space dimension due to Sobolev embeddings. The focus is on the final data case (3) there. In Section 4 we show numerical tests with the discussed fixed point iteration. We provide reconstructions of the pair ( f 1 , f 2 ) for known φ i , r i , w, as well as of the pair (φ 1 , φ 2 ) for known f i , r i , w in (1). Here, both settings of final time data (3) and of time trace data (4) are considered.

Preamble
We start with a short revision of the results in [16], since the present paper builds on these. As a matter of fact, in [16], we considered a scalar problem with only one unknown function f but the setting there is slighty more general than the one described in the introduction (and most of the rest of this paper) in the sense that a subdiffusion equation is considered. More precisely, in [16] we treat the semilinear (sub)diffusion initial boundary value problem D α t u(x,t) + Lu(x,t) = f (u(x,t)) + r(x,t) (x,t) ∈ Ω × (0, T ) ∂ u ∂ ν + γu = 0 on ∂ Ω × (0, T ) where D α t denotes the Djrbashian-Caputo fractional time derivative of order α ∈ (0, 1) which is defined by In case α = 1, D α t is the usual first time derivative and (5) becomes a semilinear parabolic equation. For details on fractional differentiation and subdiffusion equations, we refer to, e.g., [1,2,3,20,28,29], see also the tutorial on inverse problems for anomalous diffusion processes [12]; For well-posedness of the forward problem (5) with Lipschitz continuous f we refer to [11,Section 3], see also [16,Section 2].
In [16], the inverse problem of recovering the nonlinearity f from final time data g(x) ≡ u(x, T ), x ∈ Ω (assuming that all other coefficients in (5) are known) was approached by a fixed point scheme for the operator T defined by the identity (T f )(g(x)) = D α t u(x, T ; f ) − Lg(x) − r(x, T ) x ∈ Ω . where for given f , the function u(x,t; f ) denotes the solution to (5). This defines T f just on the range of g and therefore it is crucial to assume that all values of u that get inserted into f belong to this range J := g(Ω) ⊇ u(Ω × [0, T )) . (8) More precisely, it suffices to impose this condition at the solution ( f ex , u ex ) of the inverse problem, and to project the iterates onto this range during the fixed point iteration. Obviously, these considerations can be extended to the case of partial observations g(x) ≡ u(x, T ), x ∈ ω on a subset ω ⊂ Ω, under the adapted range condition g(ω) ⊇ u(Ω ×[0, T )).
Since f is an univariate function, measurements on a higher dimensional domain Ω appear to be too much in fact and it can suffice to take observations along a curve ω connecting the points of minimal and maximal values of u(x, T ).
The following example illustrates the fact that this imposes a true constraint on the class of problems that can be expected to exhibit unique identifiability. Example 2. 1. Let ϕ be an eigenfunction of −L with corresponding eigenvalue λ > 0. Take f (u) = cu as well as u 0 = ϕ , r = 0, then u(x,t; f ) = e −(λ −c)t ϕ(x). If ϕ stays positive (which is, e.g., the case is λ is the smallest eigenvalue), and c < λ , then, with ϕ = min x∈Ω ϕ(x) ≥ 0, ϕ = max x∈Ω ϕ(x) > 0, we get for the range of u over all of Ω, and therefore the range condition will be violated.
In [16,Section 3] we proved a self-mapping property of T on a sufficiently small ball in W 1,∞ (J), as well as its contractivity in the parabolic case α = 1 for T large enough provided f is strictly monotonically decreasing, which implies exponential decay of the corresponding solution or actually its time derivative D t u, as long as r vanishes or has exponentially decaying time derivative, see also Lemma 3.1 below. Such a dissipative setting is indeed crucial for proving that the Lipschitz constant of T decreases with increasing final time T , as the following counterexample shows. Example 2. 2. Again, let ϕ(x) be an eigenfunction of −L with corresponding eigenvalue λ > 0. Take f (1) (u) = c 1 u, f (2) (u) = c 2 u as well as u 0 = ϕ , r = 0, and set u (i) (x,t) = u(x,t; f i ) which can be computed explicitely as ) − T f (2) (g(x)) = u (1) t (x, T ) − u (2) t (x, T ) = (c 1 − λ )e c 1 T − (c 2 − λ )e c 2 T e −λ T ϕ(x) .
Thus, for any combination of norms · X , · Z the contraction factor T f (1) (g)−T f (2)  The two examples above clearly show that the two aims (a) range condition and (b) dissipa-tivity (for contractivity) are conflicting, at least as long as we set r = 0 as it was done here. Thus, in order to achieve (a) we need to drive the system by means of r (or alternatively, by inhomogeneos boundary conditions). Luckily, nonvanishing r does not impact the case (b) since these inhomogenieities basically cancel out when taking differences between fixed point iterates for establishing contractivity estimates. Thus it is indeed possible to have range condition and contractivity together. However, for this it is crucial to not only use found data, but to be able to design the experiment in such a way that the data exhibits the desired properties.
The analysis in [16] is restricted to the case of a single scalar equation and to one space dimension, due to the fact that it is carried out in (Hilbert) Sobolev spaces using embeddings into W 1,∞ . The aim of the analysis here is therefore twofold. First of all, we intend to generalize the approach from [16] towards identification of nonlinearities in certain systems of reaction-diffusion equations, see Section 3. 1. Secondly, we provide an analysis of the above fixed point scheme in Schauder spaces, which allows us to work in higher space dimensions, see Section 3. 2. For the latter case we have to restrict ourselves to the parabolic setting α = 1, as carrying over the regularity theory from, e.g., [6] to the subdiffusion setting, based on results from, e.g., [5], would be a major effort in itself that goes beyond the scope of this paper. This also affects the extension of the analysis from [27,26] for time trace (4) instead of final time data (3), which is why we focus on the final time observation case in the analysis section 3.
A generalization of the scalar diffusion setting to systems of an arbitrary number N of possibly interacting states can be achieved by replacing the unknown function f in (5) by a component wise defined vector valued unknown nonlinearty f = ( f 1 , . . ., f N ), whose action is not only defined on a single state but on a possible combination of these, via known functions w = (w 1 , . . ., w N ). Possible additional known interaction and reaction terms are encapsulated in a set of (now potentially also nonlinear) functions r = (r 1 , . . . , r N ).
We thus consider systems of reaction-(sub)diffusion equations of the form D α t u i (x,t) + (L u) i (x,t) = f i (w i ( u(x,t))) + r i (x,t, u(x,t)) (x,t) ∈ Ω × (0, T ) , i ∈ {1, . . .N} , for u = (u 1 , . . ., u N ) subject to the boundary and initial conditions and Well-posedness of this forward problem with Lipschitz continuous nonlinearities f i , r i is a straightforward extension of the results from [11,Section 3], [16,Section 2].
Given the self-adjoint uniformly elliptic operator −L, e.g. − L = diag(−∇ · (A∇·), . . ., −∇ · (A∇·)) + Q(x) (12) with Q : Ω → R N×N , A ∈ R N×N symmetric (uniformly) positive definite, the functions w : I := I 1 × · · · × I N → J := J 1 × · · · × J N , r : Ω × (0, T ) × I → R N , and the data u 0 , γ i , as well as measurements on a subset ω of the domain Ω we wish to determine the unknown functions This includes both cases of identifying the reaction terms f i (by setting w i (ξ 1 , ξ 2 ) = ξ i ) and of identifying the interaction terms φ i in (1). We will abbreviate the collection of unknown functions by f = ( f 1 , . . . , f N ), noting that each individual f i might have a different domain of definition J i .
Note that this setting allows for linear and nonlinear coupling among the individual states u i via the known differential operator L (often referred to as cross-diffusion) and the known functions w i as well as r i .
Throughout the analysis we will impose the range condition cf, (8), where we assume ω to be compact to guarantee (via Weierstrass' Theorem and continuity of g i ) that I i is indeed a compact interval. Here u ex is the state part of a solution ( f ex , u ex ) of the inverse problem.
Part of the analysis is based on series expansions in terms of the eigenvalues and -functions (λ n , φ n ) n∈N of the self-adjoint elliptic operator L as well as the induced Hilbert spaceṡ H σ (Ω) = {v ∈ L 2 (Ω) : , that is equivalent to the H σ (Ω) Sobolev norm provided the coefficients of L are sufficiently regular, which we assume to be the case here.
Consider now the spatially one-dimensional setting of Ω ⊆ R 1 being an open interval (0, L), and make the invertibility and smoothness assumptions that by the Inverse Function Theorem imply that w i • g : ω → J i is bijective and its inverse is in H 2 (J i ; Ω). Thus we can define the fixed point operator T : X → X , where by with the projection P : R N → I on the compact cuboid I, i.e., P i ξ = max{g i , min{g i , ξ }}.
Besides the resonstruction problem (9), (13) with final time data, that will be discussed in detail in the convergence section 3, in the numercial reconstruction section we will also consider an analogous inverse problem of recovering f = ( f 1 , . . ., for some x 0 ∈ ∂ Ω, cf. Pilant and Rundell [27,26] for the scalar case. Under the invertibility condition we can, analogously to [27,26], define a fixed point operator T : The crucial estimates of (L u) i (x 0 ,t; f ) C 0,1 (Ω) required for establishing self-mapping and contraction properties of T on a ball in X as in [27,26] could in principle like there be based on the implicit representation of u by means of the Green's function K and the solution ψ(x,t) = u(x,t; 0) of the linear problem obtained by setting f ≡ 0, and replacing r(x,t, u(x,t)) by r(x,t, 0), together with regularity estimates on the Green's function K. For Green's functions for systems and their regularity in the parabolic case α = 1, see, e.g. [6, Theorem 1, Chapter 9]. In the subdiffusion case α < 1 the Green's function is defined by the Fox H-functions, cf., e.g., [5].

Convergence of a fixed point scheme for final time data
We will now consider convergence of the fixed point scheme defined by the operator T defined by (17) for reconstructing the reaction and interaction functions f i in (9). To some extent we can here build on previous work for the case of one scalar equation and a single function f to be reconstructed in [16]. However, it is also clear that the interaction among several states can complicate the situation considerably and lead to phenomena that would not be possible with single uncoupled equations. Our aim in Section 3.1 is to explore conditions for a scenario that would allow to make some statements on self-mapping and contractivity of T. We are aware of the fact that this is far from capturing the whole multitude of possibilities and interesting cases that can arise in systems. The theoretical results are illustrated by some examples of 2 × 2 systems arising in systems biology.
Another extension made in this section is to get rid of the restriction on the spatially one dimensional setting that had to be imposed in [16] in order to enable certain Sobolev embed-dings, in particular at the transition from H s (Ω) for the space dependent function f (g(x)) to W 1,∞ (J) for the univariate function f (u). We do so by carrying out the analysis in Schauder spaces instead, which basically allows to use the same differentiability order for f (g(x)) and f (u), independently of the dimension of Ω. Since the required regularity results on PDE solutions are so far only available in the literature for α = 1, we restrict ourselves to this case in Section 3. 2. 3.1 Vectorization of the results from [16] in the spatially one-dimensional case Analogously to the proof of [16, Theorem 3.1] we can establish T as a weakly * continuous self-mapping on a sufficiently large ball in X . (15) hold and assume that κ as well asρ := sup ζ ∈I D t r(·, ·, ζ ) L Q * (0,T ;L 2 (Ω)) are sufficiently small. Then for large enough ρ > 0 the operator T defined by (17) is a self-mapping on the bounded, closed and convex set and T is weakly* continuous in X as defined in (16). Thus T has a fixed point in B.
Moreover, in the parabolic case α = 1, contractivity of T for sufficiently large final time T follows as in [16,Theorems 3.2,3.3] from the fact that for f (1) The factor D t u (2) appearing in the right hand side of this PDE decays exponentially under certain conditions. Lemma 3. 1. Let L be of the form (12) where for some c Q > −λ , with λ > 0 the smallest eigenvalue of −∇ · (A∇·) with homogeneous impedance boundary conditions, (note that here the definition of nonnegative definiteness of a matrix here does not necessarily include its symmetry) where and assume that there exist constants C r , c r > 0 such that |D t r(x,t, ξ )| ≤ C r e −c r t x ∈ Ω , t > 0 , ξ ∈ I , and that D t u i (0) = (L u 0 ) + f i (w i ( u 0 )) + r i (0, u 0 ) ∈ L ∞ (Ω).
We mention in passing that for the proof of Lemma 3.1 it was essential to work with homogeneous boundary conditions.
Note that this result does not rely on any Sobolev embeddings and therefore remains valid for higher space dimensions, provided we can make sense of the condition (26), which in one space dimension easily follows even in the general form from uniform strict monotonicity of the functions g i and x → u (2) i (x,t) for all t ≥ 0 as well as the range condition (14).

Some Examples
We now provide some examples of typical systems to which the analysis above applies.
where w i (ξ 1 , . . . , ξ N ) = ξ i , and (15) is satisfied for strictly monotone and H 2 smooth g i , with J i = I i = g i (ω). The Jacobian of w is just the identity matrix and M(x,t, (22) is satisfied for monotonically decreasing functions f j , or if the f ′ j have arbitrary sign but their positive values on I j are dominated by λ min (Q(x)) + λ . This clearly extends to , with sufficiently small interaction terms r i . A particular case of interest here is the so-called competing species interaction where r i ( u) = −u i ∑ j =i β i j u j with nonnegative coefficients β i j , cf., e.g., [13] and see Example 3.2 below for the case N = 2.
We now turn to some examples of 2 × 2 systems, where condition (22) can be verified by applying the simple criterion , and we expect both u and v to be nonnegative, thus set for some upper bounds u max , v max > 0. Nonnegativity of −M via (28) is equivalent to A sufficient condition for this to hold is i.e., a restriction on the size of the interaction β u.v, in terms of the reactions see, e.g., [17, equations (1), (2), page 61], where typically r 2 (u, v) = −r 1 (u, v) or at least the individual corresponding terms have opposite sign. More precisely, they are often of the form with nonnegative constants k i± , ν i j , µ i j , which would be the characteristic for, e.g., massaction kinetics. In view of the practically relevant setting of nonnegative states, we will again consider the nonlinear functions on In this example, , so that nonnegativity of −M via (28) is equivalent to In case of the particular form (32), condition (33) is equivalent (by considering the asymptotics as u, v → 0) to ν 21 ≥ µ 21 and ν 22 ≥ µ 22 and ν 22 k 2+ u The requirements of f ′ ∈ L ∞ ([0, u max ]), and r 1, lead to further restrictions on the exponents ν i j , µ i j to avoid singularities at vanishing u, v. Example 3. 4.
together with (42) does not impose additional conditions on the exponents in (32).

Analysis in Schauder spaces in the parabolic case and higher space dimensions
In the parabolic case α = 1, the availability of regularity results in Schauder spaces C k,β (Ω× (0, T )), cf., e.g, [6], allows to work in higher space dimensions Ω ⊆ R d , d > 1. Thus these results are new as compared to those in [16] even in case of a single PDE. Note that the Schauder space setting has already been used in [27,26] for the same nonlinearity identification problems in case of time trace (instead of final time) observations.
For simplicity of exposition we here consider the scalar case (abbreviating D t u by u t since there will be no further subscripts here) of recovering f : Note that analogously to Section 3.1, this can be extended to the system setting (9)-(13), provided a higher dimensional (with respect to space) version of the condition (15) holds and guarantees invertibility as well as smoothness of the mapping w i • g : ω → J i .

Self-mapping fixed-point operator on spaces of Lipschitz continuous functions
In order to work in function spaces over a fixed interval J := g(ω) = [g, g], like in [16] we project the values of u onto J, which can as well be written by means of a superposition operator (Pu)(x,t) = max{g, min{g, u(x,t)}} = Φ(u(x,t)) . Note that Φ is contained in W 1,∞ (R) = C 0,1 (R) which will be sufficient for the proof of T being a self-mapping on X = C 0,1 (J). Later on, when using higher order Schauder spaces to show contractivity of T, the lack of additional smoothness of Φ will remain an issue.
Moreover, we assume the range condition to hold for any exact solution ( f ex , u ex ) of the inverse problem (43), (44).
Thus we define the fixed point operator T : X → X by where for some j ∈ C 0,1 (J), the function u(x,t; j) solves Equation (46) indeed uniquely determines f + := T f on J if, e.g, ω is curve in Ω along which g is strictly monotone. Otherwise, the transition from the multivariate function f + (g) : ω → R, defined on the d-dimensional domain ω , to the real function f : J → R, defined on an interval J ⊆ R, can be carried out by metric projection, cf. [16, Lemma 3.1, Remark 3.1]. We therefore redefine T : X → X by and assume that the mapping P g : C 0,1 (ω) → C 0,1 (J) satisfies the compatibility condition for an exact solution ( f ex , u ex ) of the inverse problem (43), (44), and is continuous as a mapping from C 0,1 (ω) to C 0,1 (J), i.e., for all y ∈ C 0,1 (ω), the bound holds. If, e.g, ω is a curve with a regular parametrization ω = {x(τ) : τ ∈ [0, 1]} such that ∇g(x(τ)) ·ẋ(τ) ≥ 1 C(g) |ẋ(τ)| for all τ ∈ [0, 1], then for P g simply defined by (P g w)(g(x)) = w(x), x ∈ ω , the estimate (50) follows from Assumptions (45) and (49) imply that the f part of any such solution ( f ex , u ex ) is a fixed point of T.

Contractivity in higher order Schauder spaces
In order to prove contractivity of T, we need to move on to higher order Schauder spaces X ⊆ C 2,β (R) for some fixed β ∈ (0, 1]. The PDE estimates we will use for this purpose will rely on the Schauder space regularity theory for parabolic equations from [6]. Note that in view of the counterexample [7,Problem 4.9], the Hölder exponent β needs to be strictly positive. A first attempt to circumvent the lack of higher smoothness of the projection operator would be to replace Φ in (47) by a smoothed version Φ ε . However, in order to prove convergence of the resulting approximation as ε → 0, we would need uniform boundedness of Φ ′ ε in C 0,β (R), which -as can be readily checked -is not possible, though. Note that this lack of smoothness of Φ was not an issue in [16], where we worked in Sobolev spaces, since the superposition operator induced by Φ ′ is Lipschitz as an operator from L q to L p for any 1 ≤ p < q < ∞, cf. [9].
We therefore achieve confinement to the observable range J = g(ω) in a different manner, namely by definition of the spaces in which we work as (52) Thus f ex ∈ X implies the necessary conditions f ex ′ (g) = f ex ′ (g) = f ex ′′ (g) = f ex ′′ (g) = 0 on the exact solution f ex . To circumvent these conditions, one might, in place of just enforcing constant values outside J, impose C 2,β smooth extrapolation by defining X = { j ∈ C 2,β (R) : j| [g,∞) ∈ Π([g, ∞)) , j| (−∞,g] ∈ Π((−∞, g])} for some low dimensional spaces Π([g, ∞)), Π((−∞, g]) of polynomials or rational functions. However, this would not allow to estimate the global C 2,β (R) norm of f by its corresponding norm on the observable part J. We therefore remain with the space X defined by (52). Note that for all j ∈ X the identity j C k,β (J) = j C k,β (R) k ∈ {0, 1, 2} holds.
In order to map from final time data defined on ω to functions defined on J (or actually on all of R here), we will again use an operator P g : C 2,β (ω) → X and in place of (50), assume bounds with respect to the stronger spaces to hold, i.e., for all y ∈ C 2,β (ω), Summarizing, we define T : X → X by (48) with u(x,t, j) denoting the solution of differently from the previous subsection, where we used (47).
Contractivity for small u ex,t : We will now prove that T is a self-mapping on a C 2,β ball with appropriately chosen radius ρ > 0, and that it contracts the error f − f ex C 2,β (J) . To this end, we use the definition (48) together with (53) and the fact that the function defined by z(x,t) := u t (x,t; f ) − u ex,t (x,t) on Θ := Ω × (0, T ) satisfies the parabolic initial boundary value problem Here u = u(·, ·; f ) solves (54) with j = f , andû = u − u ex solveŝ Using (53) as well as the identity (51) we obtain that To estimate z(T ) C 2,β (ω) , we apply [6, Theorem 6, page 65], together with the fact that for the space-time cylinder Θ = (0, T ) × Ω, the space C 0,β (Θ) with the norm holds for all j ∈ C 0,β (R), w ∈ C 0,β (M) and either But the latter estimate does not yield contractivity, since for this we would need an estimate by a small multiple of f − f ex C 2,β (J) .
Thus we split z = z r + z 0 into a part z r satisfying the inhomogeneous PDE with homogeneous initial conditions z r (x, 0) = 0 and a part z 0 satisfying the homogeneous To estimate the z r , we apply [6, Theorem 6, page 65] which yields where U Z is as above and Here hence, taking the supremum over τ ∈ [0,t] on both sides i.e., Likewise, for σ 0 > d/2 + β , we have where we can avoid the potential singularity of Ψ(t; σ 0 , λ 1 ) at t = 0 in (58) by assuming which still admits the three-dimensional space case. (Note that this avoidance is not possible for Ψ(t; σ , λ 1 ) due to the requirement σ > d/2 + 2 + β made above.) However, a problem occurs in (61), since due to the singularity at τ = 0 of Ψ(τ; σ , λ 1 ), the factor sup τ∈[0,t] Ψ(τ; σ , λ 1 ) is not finite. Recall that the time dependence of f ex ′ (u ex ) led to the convolution part z 02 that forced us to take the supremum over τ ∈ [0,t] to arrive at (60). Thus, the part z 0 of z corresponding to the initial condition cannot be controlled in this setting and we need to remove it by assuming the initial condition ( f − f ex )(u 0 ) to vanish. This is in line with existing results on decay of solutions to autonomous equations, see, e.g., [19], that require time periodicity of the coefficient -in our case f ′ ex (u ex ), which is not available here, though.
Doing so, we end up with an estimate of the form provided f − f ex C 2,β (J) ≤ ρ small enough. This yields self-mapping and contractivity on a ball of radius ρ around f ex in if u ex,t C 0,β (Θ) is sufficiently small. Theorem 3. 4. Let 2 > d/2 + β , and (53) hold, and assume that f ex ∈ X defined as in (64) Then there exist ρ , κ > 0, such that for u ex,t C 0,β (Θ) ≤ κ , the operator T is a self-mapping hold for some q ∈ (0, 1) and any f 0 ∈ B X ρ ( f ex ).
Note that assuming Note that this result is valid in space dimensions d ∈ [1, 4 − 2β ), so in particular also for d = 3, as long as β < 1 2 .
Contractivity for monotonically decreasing f : To avoid the smallness assumption on u ex,t C(Θ) in Theorem 3.4, we can make use of exponential decay of u t in case of monotonically decreasing f and exponentially decaying r t .
To this end, we assume existence of a nonnegative, only space dependent potentialq, that is sufficiently close to the space and time dependent function f ex and rewrite (55) as z t − Lz +qz = ( f ex ′ (u ex ) +q)z + y u t =: where Moreover, we take advantage of the fact that without the initial data term we gain regularity of f − f ex by iterating, i.e., by (57), we have which puts us into the favorable position of only having to estimate the C 1,β norm instead of the C 2,β norm of (T f − f ex )(g) = z(T ). Note that the choice (63) that enabled us to work in higher space dimensions required us to use the H 2 norm of ( f − f ex )(u 0 ), which we estimated by the C 2 norm of f − f ex . To avoid the term f − f ex C 2 (J) in the right hand side of (61), we therefore impose ( f − f ex )(u 0 ) = 0, which can, e.g., be achieved by assuming f (0) = f ex (0) and u 0 ≡ 0.
To achieve contractivity for large enough final time T , we will additionally assume that L is of the form that r t is exponentially decaying, and that f is monotonically decreasing f ′ ≤ 0, which altogether implies that u t appearing in the right hand side of (55) decays exponentially as in [16,Section 3.3] and Lemma 3.1 above.

Thus withč 1 sufficiently small so that
we get hence, via (53) and (68), Lipschitz continuity of T 2 with a factor that decays exponentially with T .
Here it is important to note that K in (56), (68) can be chosen as independent of T in the dissipative setting we are considering here, due to the following lemma, whose prove can be found in the appendix. for someq ∈ C 2 (Ω) (depending on x only),q ≥ 0, for some µ <λ 1 , p ∈ [1, ∞), whereλ 1 > 0 is the smallest eigenvalue of −L +q.
Then there exists a constantK independent of T such that for any j ∈ C 0,β (Θ) the solution z of with homogeneous initial and conditions satisfies In particular, the assumptions of this Lemma are satisfied if c = − f ′ ex (u ex ) ≥ 0.

Reconstructions
We will show the results of numerical experiments using the basic versions of the iterative schemes defined by (17) and (21) for each of the two data types: time trace data consisting of the value of h(t) := u(x 0 ,t) for t ∈ [0, T ]; final time data g(x) := u(x, T ) for some chosen value of T . The numerical results presented will be set in one space dimension although there is no limitation in this regard (other than computational complexity of the direct solvers) as our unknowns are functions of a single variable. Also, in this setting the graphical illustrations are more transparent. Note that in one dimension the curve ω ⊂ Ω becomes Ω itself which we take to be the unit interval. We will also consider only two equations as this case encompasses most of the features of a larger system. For notational convenience we use u and v for the dependent variables in the two equations in the system. As data we took two differing initial values u 0 (x) and v 0 (x) and as boundary conditions we used (homogeneous) Dirichlet at the left endpoint and Neumann at the right; typically different for each of u and v.
We outline below the main steps used to compile the reconstruction examples shown throughout this section. 1. The domain and solvers used. (c) For the spatial case of g meas (x) this filtering used an eigenfunction basis of the elliptic operator to take into account the boundary conditions and projecting onto a basis set of its eigenfunctions using H 2 smoothing of its coefficients to obtain g(x). In the temporal case of h meas (t) where the only constraint is with the initial data function at t = 0, either a H 1 Tikhonov penalty term or a smoothing spline routine was used.
(d) Note that the iteration schemes below themselves contains no specific regularization (although the projection of each iteration onto the range of the data could be considered in this light). For more details on data smoothing and propagation of the noise through the fixed point iteration, we refer to [16, Section 3.5]. 3. Algorithm for reconstructing f 1 , f 2 in section 4.1: (a) Set f 0 1 , f 0 2 to some initial guess then for k = 0, 1, 2, . . .

Reconstructions of f 1 and f 2
In this first group of reconstructions we seek the recovery of the reaction terms f 1 (u) and f 2 (v) and assume the interaction terms between them are just given by a multiple of φ i (w) = w = uv. Our equations are then representing a "competing species" model if β < 0 and a "symbiotic relationship" if β > 0.
The magniture of β represents the strength of the coupling. The source terms r u (x,t, u) and r v (x,t, v) are assumed known if present. We used a homogeneous Dirichlet boundary conditions at x = 0 and Neumann conditions forcing in flux at x = 1. The initial conditions u 0 and v 0 were different as were r u and r v , The following sample functions to be reconstructed were used: The first of these is a version of the Zeldovich combustion model with chosen parameters that are physically relevant, the second is chosen to offer more challenge to the reconstruction process. Depending on the driving boundary conditions and the strengths of the interaction terms φ i the range over which one must recover f i can be considerable. This offers challenges from a computational viewpoint and of course if one "knew" that the correct answer was a polynomial function all the reconstruction process would be nothing other than a least squares fit in some appropriate norm to obtain a small number of constants. We are looking beyond this here and hoping to be able to detect features in these reaction and coupling terms that might drive the model rather than be purely derivative from it.
The iteration schemes (17) and (21) can be implemented pointwise or by representing the unknowns in a set of basis functions. The choice of the latter is important. While many standard models use only low degree polynomials to represent the modelling of f this is clearly a severe limitation. Using high degree polynomials is out of the question due to the severe ill-conditioning recovering Taylor coefficients from data far from the initial point. Rational functions may seem to be a good choice but again their range of accuracy is limited when used over a wide interval. In addition, this is a nonlinear fitting problem -it is also unstable under extrapolation as it is again analytic continuation. There are also other negative effects. Our initial guesses for f i may be quite distant from the actual and in this situation during the iteration process it was frequently found that the denominator of the rational function has zeros sufficiently near to the real axis that reconstructed functions f k of large amplitude resulted at a local point. This had the effect of causing failure of the iteration scheme to converge. Indeed basis functions designed for a fixed and relatively narrow range tend to be suboptimal in this setting.
We found a good choice to be moving Gaussian basis functions. These are known to effectively model many non-linear relationships and since each basis function is non-zero over a small interval this localization is useful in the current situation where f 1 (u) and f 2 (v) are totally locally defined. As a secondary consideration here this locality property results in a sparse matrix that leads to much faster computation of this phase.
As noted earlier, we are initially taking the interaction functions {φ 1 (w), φ 2 (w)} to be a constant β times the identity and also w = uv. If β = 0 then this case is just a complete decoupling of the system and the results of [16] show a unique recovery through the resulting contraction mappings. For β sufficiently small, the analysis of section 3 then shows the same result and the question becomes if this holds true for all β . Our analysis does not cover this case and as we will see below this answer certainly appears to be negative and is illustrated graphically in Figure 1 below.
The norms shown in these figures are discrete L 2 norms of the functions f at the 100 stored values as described previously. We also show reconstructions of f 1 and f 2 achieved after a given number of iterations in Figures 2 and 3.  In these figures we show the exact function as a black dashed line and the iterations in bold lines with the ordering: yellow, orange, red, light green and dark green. Typically          these correspond to the second, fourth, sixth, eighth and tenth but in he case of the slower converging scheme with β = 1.3 these are at 1, 3, 6, 9, 12. Note for the larger β values the reconstruction progresses by improvements of the values of smaller magnitude first due to the causality inherent in the time trace situation as opposed to one from final time measurements.
We should expect in fact that the convergence rate will decrease with increasing β and this is entirely borne out by numerical computations As expected, there is also a difference depending on the sign of β : the case of β < 0 corresponds to a negative interaction and might be expected to be a more stable case and turns out to be the situation. The case β > 0 indicates symbiosis between the quantities represented by u and v and should be expected to be a more unstable situation. This is certainly the outcome in the ordinary differential equation version and sufficiently strong reaction terms will dominate the diffusion effect supplied by the elliptic operator. Thus one of the two of u or v should dominate and lead to blow up of the solution -perhaps in finite time. The conditions for this are quite complex. In the case of an ordinary differential equation the nonlinear right hand side terms must,          in general, be uniformly Lipschitz and if they are positive and have asymptotic growth of order u δ for δ > 1 then this occurs. Under the effects of a diffusion operator this is more complex; the power law growth is coupled to, among other things, the dimension of the space, see [10,24,32]. The limitation we must impose in this case is that the measurements taken be up to or at a time T before the onset of the solution failing to exist. Of course, long before this situation the iteration scheme would most likely have failed. This is apparent from Figure 1  However, and we stress this part, we make no prior assumptions on any of the functions f i or φ i of analyticity or indeed any form that would be non-local in nature; the possibility that these reaction terms could have very different properties in different ranges of their arguments is a central hypothesis. Thus simply assuming that small time measurements leading to small values of u and v and the associated reconstruction of f i and φ i over this range extends to a global reconstruction is too restrictive in many physical situations.
Returning to Figure 1 the topmost pair shows the decrease in f (n) i − f i for i = 1, 2 using time trace data on the right boundary point where the solution is being driven. For both f 1 and f 2 the initial approximation was the zero function showing that a good initial approximation isn't critical here. For β greater than about 1.3 the scheme no longer converged. The bottom pair shows the analogous result for final time data. The scheme converged for β less than about 0.3 with the time measurement taken to be T = 1, but if this were reduced, somewhat larger values of β can be used. The case of β = 0.5 and T = 0.75 which approximately corresponds to the maximum value for T with is β , illustrates this situation.
The actual final reconstructions obtained from the above choices of β differ slightly. But as Figure 1 shows the individual iterations show considerable differences. Note also that the range of values of u and v depend quite strongly on β . The reconstructions obtained are shown in Figure 2 for the case of time trace data and in Figure 3 for the case of final time data with parameters as described above.

Reconstructions under added noise
We show below similar reconstructions but under noise in the data. Although from a technical perspective of considering only the norms of the unknowns and the data space these recovery problems appear to be only mildly ill-conditioned this does not take into account the relevant constants that are strongly influenced by the nonlinearities in the model. Thus while good reconstructions can be obtained with a few percent added uniform noise, much more leads to quite poor reconstructions. In Figure 4 we show the same functions f 1 (u) and f 2 (v) reconstructed under noise levels of 0.1% and 1%. These show the quite obvious distinction between the previous reconstructions where a very mild filter was added of the type noted above to offset any inaccuracy in the direct data due to truncation error in the direct solver and filtering to remove active noise.

Reconstructions of the interaction terms φ i
In this section we assume that both f 1 and f 2 are known, but φ 1 (w) and φ 2 (w) in are unknown and have to be recovered from either a pair of time trace or a pair of final time data at t = T . The sample values to be reconstructed are       where β u and β v are constants that will be used to test how large these functions can be and the iteration scheme still converge. The known functions in the model were set to We show the reconstructions of φ 1 and φ 2 graphically by again displaying four pairs of reconstruction iterates, for each of the values β = −1, 0.1, 1, 10, in Figure 5. Of course different values of β give different solution with differing ranges and this is clearly shown. The exact value is the dotted curve, the curve in yellow is the first iterate, the one in orange is the second, and the one in red the third or fifth. This shows the rapid convergence for a wide range of β values.

Epilogue
It is tempting for authors to wonder how a paper will be received and in this case the answer to the question is likely to depend on the community to which the reader belongs.
The practitioners might feel not enough attention was given for complete answers to specific problems or the range of problems was insufficient; "why was . . . not tackled? Mathematicians might have liked to see theorems containing "sufficently small" conditions replaced by estimates with tangible values. The inverse problems community seeing the complexities arising from the unknown ranges and nonlinearities themselves, might reflect, "but linear inverse problems/equations behave even better." In some sense these are valid statements. However, reaction diffusion systems are able to model an enormous range of physical problems and coupled systems of nonlinear equations are always going to impose mathematical difficulties. There are indeed easier inverse problems, but it is the above ubiquity and challenges that make them compelling.
We hope to continue this work by expanding the range of questions posed, by looking for better analytic tools and superior computational methods. There is much, much more still to be said.