Biological invasions: Deriving the regions at risk from partial measurements

Mathematical Biosciences 215 (2008) 158–166

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Biological invasions: Deriving the regions at risk from partial measurements Michel Cristofol a,b, Lionel Roques c,* a b c

LATP, CNRS/UMR 6632, CMI, Université de Provence, 39 rue Joliot Curie, 13453 Marseille Cedex 13, France Université d’Aix-Marseille III, IUT de St Jérôme, France INRA, UR546 Biostatistique et Processus Spatiaux, Domaine St Paul Site Agroparc, F-84914 Avignon, France

a r t i c l e

i n f o

Article history: Received 16 November 2007 Received in revised form 18 June 2008 Accepted 11 July 2008 Available online 24 July 2008 Keywords: Reaction–diffusion Biological invasions Inverse problem Habitat configuration Carleman estimates Simulated annealing

a b s t r a c t We consider the problem of forecasting the regions at higher risk for newly introduced invasive species. Favourable and unfavourable regions may indeed not be known a priori, especially for exotic species whose hosts in native range and newly-colonised areas can be different. Assuming that the species is modelled by a logistic-like reaction–diffusion equation, we prove that the spatial arrangement of the favourable and unfavourable regions can theoretically be determined using only partial measurements of the population density: (1) a local ‘spatio-temporal’ measurement, during a short time period and, (2) a ‘spatial’ measurement in the whole region susceptible to colonisation. We then present a stochastic algorithm which is proved analytically, and then on several numerical examples, to be effective in deriving these regions. Ó 2008 Elsevier Inc. All rights reserved.

1. Introduction Because of trade globalisation, a substantial increase in biological invasions has been observed over the last decades (e.g. Liebhold et al. [1]). These invasive species are, by definition [2], likely to cause economic or environmental harm or harm to human health. Thus, it is a major concern to forecast, at the beginning of an invasion, the areas which will be more or less infested by the species. Because of their exotic nature, invading species generally face little competition or predation. They are therefore well adapted to modelling via single-species models. Reaction–diffusion models have proved themselves to give good qualitative results regarding biological invasions (see the pioneering paper of Skellam [3], and the books [4,5] and [6] for review). The most widely used single-species reaction–diffusion model, in homogeneous environments, is probably the Fisher–Kolmogorov [7,8] model:

ut ¼ DDu þ uðl  cuÞ;

t > 0; x 2 X  RN ;

ð1:1Þ

where u ¼ uðt; xÞ is the population density at time t and space position x, D is the diffusion coefficient, l corresponds to the constant intrinsic growth rate, and lc is the environment’s carrying capacity. Thus c measures the susceptibility to crowding effects. * Corresponding author. Fax: +33 432722182. E-mail addresses: [email protected] (M. Cristofol), lionel.roques@ avignon.inra.fr (L. Roques). 0025-5564/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2008.07.004

On the other hand, the environment is generally far from being homogeneous. The spreading speed of the invasion, as well as the final equilibrium attained by the population are in fact often highly dependent on these heterogeneities [4,9–11]. A natural extension of (1.1) to heterogeneous environments has been introduced by Shigesada, Kawasaki, Teramoto [12]:

ut ¼ rðDðxÞruÞ þ uðlðxÞ  cðxÞuÞ;

t > 0; x 2 X  RN :

ð1:2Þ

In this case, the diffusivity matrix DðxÞ, and the coefficients lðxÞ and cðxÞ depend on the space variable x, and can therefore include some effects of environmental heterogeneity. In this paper, we consider the simpler case where DðxÞ is assumed to be constant and isotropic and c is also assumed to be positive and constant:

ut ¼ DDu þ uðlðxÞ  cuÞ;

t > 0; x 2 X  RN :

ð1:3Þ

The regions where l is high correspond to favourable regions (high intrinsic growth rate and high environment carrying capacity), whereas the regions with low values of l are less favourable, or even unfavourable when l < 0. In what follows, in order to obtain clearer biological interpretations of our results, we say that l is a ‘habitat configuration’. With this type of model, many qualitative results have been established, especially regarding the influence of spatial heterogeneities of the environment on population persistence, and on the value of the equilibrium population density [4,9,13–15]. However, for a newly introduced species, like an invasive species at the beginning of its introduction, the regions where l is high or low

159

M. Cristofol, L. Roques / Mathematical Biosciences 215 (2008) 158–166

may not be known a priori, particularly when the environment is very different from that of the species native range. In this paper, we propose a method of deriving the habitat configuration l, basing ourselves only on partial measurements of the population density at the beginning of the invasion process. In Section 2, we begin by giving a precise mathematical formulation of our estimation problem. We then describe our main mathematical results, and we link them with ecological interpretations. These theoretical results form the basis of an algorithm that we propose, in Section 3, for recovering the habitat configuration l. In Section 4, we provide numerical examples illustrating our results. These results are further discussed in Section 5.

~ be a function in M, and let v ~ be the solution of the linear parLet l abolic problem ðPl~ ;0 Þ. We define a functional Gl , over Rþ  M, by

~ Þ ¼ kot uc  ot v ~k2L2 ððt ;t ÞxÞ þ kDuc ðT 0 ; Þ  Dv ~ðT 0 ; Þk2L2 ðXÞ Gl ðc; l 0 1 ~ðT 0 ; Þk2L2 ðXÞ ; þ kuc ðT 0 ; Þ  v where uc is the solution of ðPl;c Þ. This functional Gl quantifies the ~ on the set where uc has been measured. gap between uc and v Theorem 2.2. The functions

~ k2L2 ðX Þ 6 kl  l 1

2. Formulation of the problem and main results

l; l~ 2 M being given, we have:

C ~ Þ; Gl ð0; l ui 2

~ 2 M and for some positive constant C ¼ CðX; X1 ; for all l x; Be ; D; t0 ; t1 ; ui =ui Þ.

2.1. Model and hypotheses We assume that the population density uc is governed by the following parabolic equation:

8 > < ot uc ¼ DDuc þ uc ðlðxÞ  cuc Þ; uc ðt; xÞ ¼ 0; t > 0; x 2 oX; > : uc ð0; xÞ ¼ ui ðxÞ in X;

2.3. Estimating the habitat configuration

t > 0; x 2 X; ðPl;c Þ

where X is a bounded subdomain of Rd with boundary oX. We will denote Q :¼ ð0; þ1Þ  X and R :¼ ð0; þ1Þ  oX. The growth rate function l is a priori assumed to be bounded, and to take a known constant value outside a fixed compact subset X1 of X:

l 2 M :¼ fq 2 L1 ðXÞ; M 6 q 6 M a:e:; and q  m in X n X1 g; for some constants m; M, with M > 0; the notation ‘a.e.’ means ‘almost everywhere’, which is equivalent to ‘except on a set of zero measure’. The initial population density ui ðxÞ is assumed to be bounded (in C 2 ðXÞ), and bounded from below by a fixed positive constant in a fixed closed ball Be  X1 , of small radius e:

D :¼ f/ P 0; / 2 C 2 ðXÞ; k/kC2 ðXÞ 6 ui ; / P ui in Be ; g;

ð2:4Þ

for some positive constants ui and ui . Absorbing (Dirichlet) boundary conditions are assumed. Remark 2.1. Absorbing boundary conditions mean that the individuals crossing the boundary immediately die. Such conditions can be ecologically relevant in numerous situations. For instance for many plant species, seacoasts are lethal and thus constitute this kind of boundaries. For technical reasons we have to introduce the subset X1 , such that, in the interface between X1 and X, l takes a known value m. This value is typically negative, indicating that, near the lethal boundary, the environment is unfavourable. This assumption is not very restrictive since, in fact, X1 can be chosen as close as we want to X. For precise definitions of the functional spaces L2 , L1 and C 2 as well as the other mathematical notations used throughout this paper, the reader can refer, e.g., to [16].

The proof of this result is given in Appendix A.1. It bears on a Carleman-type estimate. Biological interpretation: This stability result means that, in the linear case corresponding to Malthusian populations (c ¼ 0), two ~ cannot lead to close populadifferent habitat configurations l; l ~. Indeed, having population densities that are tion densities u0 ; v close to each other in the two situations, even on a very small region x, during a small time period ðt 0 ; t 1 Þ, and in the whole space X at a single time T0, would lead to small Gl values, and therefore, from Theorem 2.2, to close values of the growth rate coefficients l ~. and l Theorem 2.2 implies the following uniqueness result: Corollary 2.3. If v is a solution of both ðP l;0 Þ and ðPl~ ;0 Þ, then a.e. in X1 , and therefore in X.

l ¼ l~

Biological interpretation: In the linear case (c ¼ 0), if two habitat ~ lead to identical population densities u0 , v ~, even configurations l; l on a very small region x, during a small time period ðt 0 ; t 1 Þ, and in the whole space X at a single time T0, then these habitat configurations are identical. Next we have the following result: ~ Þ  Gl ðc; l ~ Þj ¼ Oðui 3 Þ, as ui ! 0: Theorem 2.4. We have1 jGl ð0; l The proof of this result is given in Appendix A.2. Biological interpretation: Assume that the habitat configuration l is not known, but that we have measurements of the population density uc , governed by the full non-linear model (1.3). Consider a ~ in M such that the population density v ~ obtained configuration l as a solution of the linear model ðP 0;l~ Þ has values close to those ta~ Þ is close ken by the population density uc , in the sense that Gl ðc; l to 0. If the initial population density is far from the environment carrying capacity, then ui  lc, ui is small and, from Theorem 2.4, ~ Þ is also close to 0. Thus Theorem 2.2 implies that the habGl ð0; l ~ is an accurate estimate of l. In Section 3, we itat configuration l propose an algorithm to obtain explicitly such estimates of l. Remark 2.5. In fact, the term Oðui 3 Þ increases exponentially with time t1 . Thus, obtaining accurate estimates of l require, in practice, to work with small times i.e. at the beginning of the invasion. 2.4. Forecasting the fate of the invading population

2.2. Main question The main question that we presented at the end of Section 1 can now be stated: for any time-span ðt0 ; t1 Þ, and any non-empty subset x of X1 , is it possible to estimate the function lðxÞ in X, basing ourselves only on measurements of uc ðt; xÞ over ðt0 ; t1 Þ  x, and on a single measurement of uc ðt; xÞ in the whole domain X at a time 1 T 0 ¼ t0 þt ? 2

~ of l enables us to give an The knowledge of an L2 -estimate l estimate of the asymptotic behaviour of the solution uc of ðP l;c Þ, as t ! þ1, and especially to know whether the population will 1 ~ ; ui ; ui ; ui ; cÞ and gðl; l ~ ; ui ; ui ; ui ; cÞ, are written f ¼ OðgÞ as Two functions f ðl; l ~ , ui , ui , ui and c, such that g ! 0 if there exists a constant K > 0, independent of l, l j f j6 K j g j for g small enough.

160

M. Cristofol, L. Roques / Mathematical Biosciences 215 (2008) 158–166

become extinct or not. Indeed, as t ! þ1, it is known that (see e.g. [9], for a proof with another type of boundary condition) the solution uc ðt; xÞ of ðP l;c Þ converges to the unique non-negative and bounded solution pc of

DDpc ¼ pc ðlðxÞ  cpc Þ in X;

ðSl;c Þ

with pc ¼ 0 on oX. Moreover, pc  0 if and only if k1 ½l P 0, where k1 ½l is the smallest eigenvalue of the elliptic operator L : w7!  DDw  lðxÞw, with Dirichlet boundary conditions. On the other hand, if k1 ½l < 0, then pc 6 0 in X (note that c does not appear in the definition of k1 ). We have the following result. ~ n ÞnP0 in M, such that Proposition 2.6. Let us consider a sequence ðl l~ n ! l in L2 ðXÞ as n ! 1. ~c;n of the problem ðSl~ n ;c Þ converges to pc as (a) The solution p n ! 1, uniformly in X. ~ n  ! k1 ½l as n ! þ1. (b) k1 ½l The proof of this result is classical and can be found in [9,17]. Biological interpretation: Assume that the habitat configuration l is not known. We know that, in large times, the population density uc will tend to an unknown steady state pc (possibly 0, in case of extinction of the population). The part (a) of the above proposi~ of l, then tion means that, if we know an accurate (L2 -) estimate l ~c of the steady state pc , prowe can deduce an accurate estimate p vided the coefficient c is known. Part (b) shows that, even if c is not ~ of l enables to obtain an estimate of known, having an estimate l k1 ½l, and therefore to forecast whether the species will survive or not. Indeed the sign of k1 ½l controls the fate of the invading species (persistence if k1 ½l < 0 and extinction if k1 ½l P 0, see [9,13,14,17] for more details). 3. Simulated annealing algorithm Let ðt 0 ; t 1 Þ be a fixed time interval, and x  X1 be fixed. We assume that we have measurements of the solution uc ðt; xÞ of ðPl;c Þ 1 over ðt 0 ; t 1 Þ  x, and of uc ðt0 þt ; xÞ in X. However, the function l 2 and the constant c are assumed to be unknown. Our objective is to build an algorithm for recovering l. Remark 3.1. When the function uc is known, the computation of Gl ðc; Þ does not require the knowledge of c. The function l is assumed to belong to a known finite subset E of M, equipped with a neighbourhood system. We build a se^ n of N elements of E with the following simulated annealquence l ing algorithm: n¼0 ^0 Initialise l while n 6 N ^ c;n Choose randomly a neighbour m of l ^ nÞ if Gl ðc; mÞ 6 Gl ðc; l l^ nþ1 m else Choose randomly with an uniform law w 2 ð0; 1Þ Gl ðc;^ ln ÞGl ðc;mÞ

HðnÞ if w < e l^ nþ1 m else l^ nþ1 l^ n endif endif n nþ1 endwhile

The sequence HðnÞ (cooling schedule) is composed of real positive numbers, decreasing to 0. The simulated annealing ^ n of elements of E. It is known algorithm gives a sequence l (see e.g. [18]) that, for a cooling schedule HðnÞ which converges sufficiently slowly to 0, this sequence converges in ^ of Gl ðc; :Þ in E (but see Remark L2 ðXÞ to a global minimiser l 3.2). Moreover, from Theorems 2.2 and 2.4, we have

^ n k2L2 ðXÞ 6 kl  l

C C ^ nÞ 6 ^ n Þ þ e1 ðui Þ; Gl ð0; l Gl ðc; l ui 2 ui 2

where e1 is a real-valued function such that e1 ðsÞ ! 0 as s ! 0. Since l 2 E we obtain that, for n large enough:

^ n Þ 6 Gl ðc; lÞ; Gl ðc; l and, from Appendix A.1

C Gl ðc; lÞ ! 0 as ui ! 0; ui 2 the ratio ui =ui being kept constant. We finally get:

^ n k2L2 ðXÞ ! 0 as ui ! 0 and n ! þ1; kl  l for a fixed ratio ui =ui . Thus, for ui small enough, and for n large en^ n is as close as we want to l, in the L2 -sense. ough, l Remark 3.2. The cooling rate HðnÞ leading to the exact optimal configuration with probability 1 decreases very slowly (logarithmically) and cannot be used in practice; see [19] for a detailed discussion. Empirically, a good trade-off between quality of solutions and time required for computation is obtained with exponential cooling schedules of the type HðnÞ ¼ H0  an , with a < 1, first proposed by Kirkpatrick et al. [20]. Many other cooling schedules are possible, but too rapid cooling results in a system frozen into a state far from the optimal one. The starting temperature H0 should be chosen high enough to initially accept ^ nþ1 all changes l m, whatever the neighbour m. For this type of algorithm, there are no general rules for the choice of the stopping criterion (see [19]), which should be heuristically adapted to the considered optimisation problem. 4. Numerical computations In this section, in one-dimensional and two-dimensional cases, we check that the algorithm presented in Section 3 can work in practice. In each of the four following examples, we fixed the sets X, X1 and M and we defined a finite subset E  M equipped with a neighbourhood system. Then, for a fixed habitat configuration l 2 E we computed, using a second-order finite elements method, the solution uðt; xÞ of ðPl;c Þ, for D ¼ 1, c ¼ 0:1, t 2 ð0; 0:5Þ, and for a given initial population density ui . Then, we fixed t 0 ¼ 0:1, t1 ¼ 0:4, and a subset x  X1 , and we stored the values of 1 ; xÞ, for x 2 X and uðt; xÞ for ðt; xÞ 2 ðt 0 ; t 1 Þ  x. Using only uðt0 þt 2 ^ n Þ of elements of E, these values, we computed the sequence ðl defined by the simulated annealing algorithm of Section 3, the ^ 0 being sampled arbitrarily, in a uniform law, among function l the elements of E. In the following examples, we used H0 ¼ 100, and HðnÞ ¼ 100  0:99n for the cooling schedule. Our stopping criterion ^ n during 500 was to have no change in the configuration l iterations. The rigorous definitions of the sets E and of the associated neighbourhood systems which are used in the following examples can be found in Appendix A.3.

M. Cristofol, L. Roques / Mathematical Biosciences 215 (2008) 158–166

4.1. One-dimensional case Assume that X ¼ ð0; 100Þ, X1 ¼ ½10; 90, x ¼ ½55; 58, M ¼ 2, m ¼ 1 and ui ¼ 0:1ð1  x=100Þ sinðpx=25Þ2 . Example 1. the set E is composed of binary step functions taking only the values m and M. 1 The function l 2 E and the measurement uc ðt0 þt 2 ; xÞ are depicted in Fig. 1(a) and (b). Two elements of E are said to be neighbours if they differ only on an interval of length 1. ^ n Þ stabilised on the exact configuration The sequence ðl about N ¼ 1500 iterations.

l after

Example 2. the set E is composed of step functions which can take 21 different values between m and M. 1 The function l 2 E and the measurement uc ðt0 þt 2 ; xÞ are depicted in Fig. 2(a) and (b). Two elements of E are said to be neighbours if they differ by ðM  mÞ=20 on an interval of length 1.

a

161

^ n Þ stabilised on a configuration l ^ (Fig. This time, the sequence ðl 2(c)) after 7500 iterations. The mean error in this case was R 100 1 ^ ðxÞjdx ¼ 0:05: jlðxÞ  l 100 0 4.2. Two-dimensional case Assume now that X ¼ ð0; 20Þ  ð0; 20Þ, X1 ¼ ½2; 18  ½2; 18, and that x is the closed ball of centre ð7; 7Þ and radius 3. Assume also that M ¼ 2 and m ¼ 1, and that the initial data is ui ¼ 0:1xy=400 sinðx=4Þ2 sinðy=4Þ2 . Example 3. E is composed of binary functions which can only take the values m and M on each cell of a regular lattice. 1 We fixed l 2 E as in Fig. 3(a). The measurement uc ðt0 þt ; xÞ is 2 depicted in Fig. 3(b). Two elements of E are said to be neighbours if they differ only on one cell of the lattice. ^ n Þ stabilised on the exact configuration l after The sequence ðl 3000 iterations.

b

2.5 2 1.5

µ(x)

1 0.5 0 0.5 1 1.5

0

20

40

60

80

100

x Fig. 1. Example 1: (a) exact habitat configuration l; (b) measurement of uc ð0:25; xÞ which was used for recovering l; the domain x is delimited by two black dots. The exact configuration l was recovered after 1500 iterations of the algorithm.

a

b

2.5 2 1.5

µ(x)

1 0.5 0 0.5 1 1.5

0

20

40

60

80

100

x

c

2.5 2 1.5

µ(x)

1 0.5 0 0.5 1 1.5

0

20

40

60

80

100

x Fig. 2. Example 2: (a) exact habitat configuration l; (b) measurement of uc ð0:25; xÞ which was used for recovering ^ :¼ l ^ 7500 , obtained after 7500 iterations. (c) configuration l

l; the domain x is delimited by two black dots;

162

M. Cristofol, L. Roques / Mathematical Biosciences 215 (2008) 158–166

a

b

20

x2

15

10

5

0

5

10 x1

15

20

Fig. 3. Example 3: (a) exact habitat configuration l. In the black regions, the depicted function takes the value 2; in the white regions, it takes the value 1. (b) Measurement of uc ð0:25; xÞ which was used for recovering l; the domain x is delimited by a black circle. The exact configuration l was recovered after 3000 iterations.

environment carrying capacity (it can be reasonably assumed at the beginning of an invasion). In such a situation, the position of the favourable and unfavourable regions, modelled through the intrinsic growth rate coefficient l, may not be known a priori. This is especially true for exotic species whose hosts in native range and newly-colonised areas can be different. From our results, in the ‘ideal case’ considered here, the position of these regions can be obtained through partial measurements of the population density. These partial measurements consist in two samples of the population density: (1) a ‘spatio-temporal’ measurement, but very locally (in the small subset x) and during a short time period and, (2) a ‘spatial’ measurement in the whole region susceptible to colonisation (X). The stochastic algorithm presented in Section 3 shows explicitly how to reconstruct the habitat arrangement l from the above partial measurements of the population density. This algorithm was proved to be effective in both one-dimensional and two-dimensional cases, in Section 4, through several numerical experiments. In Examples 1 and 3, the algorithm converged to the exact habitat

Example 4. E is composed of functions which can take 21 different values between m and M on each cell of a regular lattice. 1 The configuration l and the measurement uc ðt0 þt ; xÞ are de2 picted in Fig. 4(a) and (b). Two elements of E are said to be neighbours if they differ by ðM  mÞ=20 on one cell of the lattice. ^ n Þ stabilised on a configuration l ^ after 14,000 The sequence ðl R 20 R 20 iterations (Fig. 4(c)). The mean error was 2012 0 0 jlðx; yÞ l^ ðx; yÞjdxdy ¼ 0:04:

5. Discussion and conclusion We have shown that, for an invasive species whose density is well modelled by a reaction–diffusion equation, the spatial arrangement of the favourable and unfavourable regions can be measured indirectly through the population density at the beginning of the invasion. More precisely, we considered a logistic-like reaction–diffusion model, and we placed ourselves under the assumption that the initial population density was far from the

a

b

20 1.5 15

x2

1 0.5

10

0 5 0.5 0

5

10 x

15

1

20

1

c

20

x2

15

10

5

0

5

10 x1

15

20

Fig. 4. Example 4: (a) exact habitat configuration l; (b) measurement of uc ð0:25; xÞ which was used for recovering l; the domain x is delimited by a black circle; (c) ^ 14000 . The darker the regions in (a), (b) and (c), the higher the values of the depicted functions. approached configuration l

163

M. Cristofol, L. Roques / Mathematical Biosciences 215 (2008) 158–166

configuration. In Examples 2 and 4, the sizes of the sets of possible habitat configurations were increased compared to Examples 1 and 3. In those cases, the algorithm converged to configurations which were close to the exact ones. It is noteworthy that the spatial measurement in X and the habitat arrangement l can have very different shapes; therefore, l cannot be straightforwardly deduced from this measurement. These results can be helpful in preventing biological invasions. Indeed, a simple protocol, consisting of placing one trap in the invaded region, and recording the number of individuals captured by this trap over a short time-period (depending on the species characteristics), and performing a single survey of the number of individuals and their position in the whole considered region should allow, from our results, to detect the favourable areas, and to treat them preventively. As we have emphasised in Proposition 2.6.a, the knowledge of an estimate of the habitat arrangement l also allows us to forecast the final population density, and therefore to detect the regions at higher risk, for instance in the case of harmful species. As recalled in Proposition 2.6.b, having a good estimate of the habitat arrangement l is also crucial to forecast the fate of the invasive species: persistence or extinction. On the other hand, we have to underline that our approach may not be adapted to some species, especially those which colonies are made of few individuals. Indeed, the diffusion operator of our model can be obtained as the macroscopic limit of uncorrelated random walks. With such an operator, by the parabolic maximum principle [16], it is known that, even with a compactly supported initial population density, the solution of our model is strictly positive everywhere on the domain as soon as t > 0. This means that the solution, and therefore the information, propagate with infinite speed, which is not realistic for discrete populations. This could induce a practical limitation of our method to a certain type of species, which are well modelled by continuous diffusion processes even at low densities (typically some insect or plant species, with high carrying capacity and growth rate). Note also that some of the mathematical tools used in this paper, and especially Carleman estimates (see Appendix A.1), were initially not adapted to the non-linear model considered here. Thus, we first considered, in Theorem 2.2, the linear – or Malthusian – case. For populations whose density is far from the environment carrying capacity, the linear and the non-linear problem have close solutions. In this situation, Theorem 2.4 extended the result of Theorem 2.2 to the non-linear case of a logistic growth. The results of this paper could be immediately extended to the case of spatially varying functions cðxÞ. Another easy extension would be, for the ‘spatio-temporal’ measurement, to use a partial boundary observation on a part Cþ of the domain boundary oX instead of sampling the population over a small domain x. Using a new Carleman estimate (see [21]) we are indeed able to write a stability result for the coefficient l, similar to that of Theorem ~Þk2L2 ððt ;t ÞCþ Þ instead of kot uc  2.2, but with kon ðot uc Þ  on ðot v 0 1 2 ~ ot vkL2 ððt0 ;t1 ÞxÞ in the definition of the functional G. Acknowledgments The authors would like to thank the anonymous reviewers for their helpful comments and insightful suggestions. This study was partly supported by the French ‘Agence Nationale de la Recherche’ within the project URTICLIM ‘Anticipation des effets du changement climatique sur l’impact écologique et sanitaire d’insectes forestiers urticants’ and by the European Union within the FP 6 Integrated Project ALARM- Assessing LArge-scale environmental Risks for biodiversity with tested Methods (GOCE-CT-2003506675).

Appendix A Let us introduce the following notations: for all t; t 0 2 R, with 0 t0 > t, we denote Q t0t ¼ ðt; t 0 Þ  X and Rtt ¼ ðt; t0 Þ  oX. Throughout this section, with a slight abuse of notation, we designate by C any upper bounds in our computations, provided they only depend on the parameters X, X1 , x, Be , D, t0 , t1 , ui =ui . A.1. Proof of Theorem 2.2 A.1.1. Carleman estimate We recall here a Carleman-type estimate with a single observation. Let b be a function in C2 ðXÞ such that

g > 0 1 < b < 2 in X; b ¼ 1 on oX; minfjrbðxÞj; x 2 X n x and on b < 0 on oX; where n denotes the outward unit normal to oX. For k > 0 and t 2 ðt0 ; t 1 Þ, we define the following weight functions

uðt; xÞ ¼

ekbðxÞ ; ðt  t 0 Þðt 1  tÞ

gðt; xÞ ¼

e2k  ekbðxÞ : ðt  t 0 Þðt1  tÞ

Let q be a solution of the parabolic problem

8 t1 > < ot q  DDq þ aðxÞq ¼ f ðt; xÞ in Q t0 ; ðPÞ q ¼ 0 on Rtt10 ; > : qð0; xÞ ¼ q0 ðxÞ in X; for some functions f 2 L2 ðQ tt10 Þ, and a; q0 2 L1 ðXÞ. Then the following results are proved in [22]: Lemma 6.1. Let q be a solution of (P). Then, there exist three positive constants k0 , C 0 and s > 1, depending only on X, x, t 0 and t1 such that, for any k P k0 , the next inequalities hold:

ðaÞ kM1 ðesg qÞk2L2 ðQ t1 Þ þ kM2 ðesg qÞk2L2 ðQ t1 Þ þ sk2 t0

3 4

þs k

Z

e

t

2sg

Q t1

þ

t0

e

t

3 4

u q 6 C0 s k

0

Z



3 2

Z

t0

#

2sg

ðf  aqÞ

2

t1

Z

t

Q t1

e2sg ujrqj2

0

Z

2sg

e

3 2

uq

x

;

Q t1 0

where M1 and M 2 are defined by M1 w ¼ DDw  s2 k2 D j rbj2 u2 wþ sðot gÞw, and M2 w ¼ ot w þ 2skDurb:rw þ 2sk2 Du j rbj2 w: Moreover,

ðbÞ s2 k2 þ s2 k2

Z

t

Z

Q t1 0

 Z e2sg u3 q2 6 C 0 s2 k2

t

Q t1 1 2

þs k

e2sg u1 ½ðot qÞ2 þ ðDqÞ2  þ

Z

0

# e

t

2sg

t1

t0

Z

Z

e2sg ujrqj2

t

Q t1 0

e2sg u3 q2

x

2

ðf  aqÞ :

Q t1 0

A.1.2. Stability estimate with one observation ~ 2 M. We consider the solutions v and e Let l; l v of the linear ~, y ¼ ot w, problems ðP l;0 Þ and ðP l~ ;0 Þ, respectively. We set w ¼ v  v ~ . The function y is a solution of: and r ¼ l  l

8 ~ > < ot y ¼ DDy þ ly þ rot v yðt; xÞ ¼ 0 > : yð0; xÞ ¼ rui ðxÞ

in Q tt10 ; on Rtt10 ; in X;

ð6:5Þ

The function gðx; tÞ attains its minimum value with respect to the 1 time at t ¼ T 0 ¼ t0 þt . We set w ¼ esg y. Using the operator M2 , intro2 duced in Lemma 6.1, we introduce, following [23] (see also [24] and [25]):

164

M. Cristofol, L. Roques / Mathematical Biosciences 215 (2008) 158–166

Z

I¼

Z

T0

t0

since wðt 0 Þ ¼ 0 and ru ¼ kurb. As a consequence, for k > 1, and since Db is bounded in X, we finally get

M 2 ww:

X

Z

Let k0 be fixed as in Lemma 6.1.

jIj 6 C s3=2 k2

Z

t1

Z

e2sg u3 y2 þ s3=2 k2

Z t

Q t1

x

t0

t0

v Þ2 : e2sg r2 ðot e

sk2

 Z e2sg uy2 6 C k2

Z 0

Q Tt

Proof. From the Hölder inequality, we have:

jIj 6 s3=2 k2

T0

t0

T0

Z

t0

X

2 2

þs k e2sg y2

1 3=2 2 k kM2 wk2L2 ðQ t1 Þ þ s3 k4 s t0 4

:

kM 2 wk2L2 ðQ t1 Þ þ 2s3 k4

t Q t1 0

t0

" 6 C s3 k 4

Z

t1

t0

Z

2

C k

0

Z

t Q t1 0

e2sg u3 y2 ;

t Q t1 0

ð6:9Þ

kM 2 wk2L2 ðQ t1 Þ þ s3 k4 t0

6 C s3 k 4

Z

t1

t0

Z

t

#

Z Z t

Q t1

x

e2sg r2 ðot e v Þ2 :

ð6:10Þ

0

The conclusion of Lemma 6.2 follows from (6.7) and (6.10).

e

2sgðT 0 ;xÞ

X

0



2

3=2 2

ðr e v ðT ; xÞÞ 6 C s

Z

t1

e

2sg

3 2

uy

Z

I¼

1 2

0 Q Tt 0

ot ðw2 Þ  skD

e2sgðT ;xÞ r2 ðot e v Þ2

0



0

e2sgðT ;xÞ ðDDwðT 0 ; xÞ þ lwðT 0 ; xÞÞ2 dx :

r  ðurbÞw2 þ 2sk2 D

0 Q Tt 0

Z 0 Q Tt 0

ujrbj2 w2 : ð6:11Þ

1 2

X

0

2

2

wðT ; Þ ¼ I  sk D

Z

2

Q Tt

0

0

e

3 2

2 2

u y þs k

#

Z

e

t

2sgðT 0 ;xÞ

Q t1

x

2

2

r ðot v~Þ :

0

and, since s > 1, the estimate of Lemma 6.3 follows.

h

~ 6 ui eMt1 , and jot v ~j 6 ðD þ MÞui eMt1 in Lemma 6.4. We have 0 6 v; v Q tt10 . Proof. From the parabolic maximum principle, we know that ~ P 0 in Q tt1 . Let h be the solution of the ordinary differential v; v 0 equation 0

h ¼ hM on Rþ ;

ð6:14Þ

The function h is increasing and Hðt; Þ :¼ hðtÞ is a supersolution of ~. As a consequence of the parabolic the equations satisfied by v and v maximum principle, we have

ð6:15Þ

in Q ; on R;

ð6:16Þ

in X:

Since ui 2 D, qð0; xÞ 2 L1 ðXÞ. Moreover, H1 ðt; xÞ :¼ ðD þ MÞui eMt and H2 ðt; xÞ :¼ ðD þ MÞui eMt are respectively sub- and supersolutions of (6.16). The parabolic maximum principle leads to the inequalities:

ðD þ MÞui eMt1 6 q 6 ðD þ MÞui eMt1 in Q tt10 :

2

u j rbj w dxdt þ skD

Z

2

Q Tt

0

0

uDbw dxdt

ð6:17Þ

h From Lemmas 6.3 and 6.4, and since follows that

Z

r vanishes outside X1 , it

0

X1

We then obtain

Z

2sg

0

t

Q t1

Proof. Using integration by parts over X and the boundary condition w ¼ 0 on Rtt10 , we get:

Z

Z

8 ~q > < ot q ¼ DDq þ l qðt; xÞ ¼ 0 > : qð0; xÞ ¼ DDui ðxÞ þ l~ ui

X

Z

t1

~. The function q satisfies: Let us set q :¼ ot v

x

t0

þ s1 k2 þ

Z

k

Z

ð6:13Þ

~ 6 hðt1 Þ ¼ ui eMt1 ; in Q tt1 : 0 6 v; v 0

h

Lemma 6.3. Let k P k0 . There exists a constant C such that

Z

2½l y þ r ðot e vÞ  :

hð0Þ ¼ ui :

Q t1

e2sg u3 y2 þ

# 2

2

If we now observe that

e2sg u3 y2 0

"

Z

t0

(

for k large enough. Combining (6.8) and (6.9), we obtain:

Z

2 2

r2 v~ðT 0 ; xÞ2 6 2yðT 0 ; xÞ2 þ 2½DDwðT 0 ; xÞ þ lwðT 0 ; xÞ2 ;

e2sg 2½l2 y2 þ r2 ðot e v Þ2  :

t Q t1 0

Furthermore, since l is bounded, and since u is bounded from below by a positive constant, independent of k, we get that

2e2sg l2 y2 6 s3 k4

e

t

e2sg u3 y2

x

we get:

#

ð6:8Þ

Z

Z

yðT 0 ; xÞ ¼ DDwðT 0 ; xÞ þ lwðT 0 ; xÞ þ rot e v ðT 0 ; xÞ;

Z

x

"

ð6:7Þ

Q t1

e2sg u3 y2

e2sg u3 y2 þ

ð6:12Þ

Arguing as in the proof of Lemma 6.2 for equation (6.9), and since, for all x 2 X, the function t7!gðt; xÞ attains its minimum over ðt 0 ; t 1 Þ at t ¼ T 0 , we finally obtain that, for k large enough, the last term in (6.12) is bounded from above by:

X

Applying inequality (a) of Lemma 6.1 to q :¼ y, we obtain that there exists C > 0, such that

Z

e2sg uy2 ;

X

0

e2sg u3 y2 :

t

2sg

Q t1

!

Z

Z

1=2

Thus using Young’s inequality, we obtain

jIj 6

t1

t0

0

1=2  Z ðM 2 wÞ2 s3 k 4

Z

0Þ for some constant C. Using u 6 ðt1 t u3 and Lemma 6.1, we get 4 that:

0

Z

T0

2

#

ð6:6Þ

Z

Z

X

Lemma 6.2. Let k P k0 . There exists a constant C such that

"

0

e2sgðT ;xÞ yðT 0 ; xÞ2 dx 6 2jIj þ Csk2

e2sgðT ;xÞ r2 ð e v ðx; T 0 Þ2  Cs1 k2 ðD þ MÞ2 ui 2 e2Mt1 Þdx

6 Cs3=2 k2

Z

t1

t0

Z x

þ wðT 0 ; xÞ2 dx:

e2sg u3 y2 þ C

Z

0

e2sgðT ;xÞ ½ðDwðT 0 ; xÞÞ2

X

ð6:18Þ

165

M. Cristofol, L. Roques / Mathematical Biosciences 215 (2008) 158–166

Lemma 6.5. The ratio ui =ui being fixed, there exists r > 0, indepen~ , ui , ui and ui , such that v ~ðT 0 ; Þ > ui r in X1 . dent of l

kot uc kL2 ðQ t1 Þ ¼ Oðkuc ðl  cuc ÞkL2 ðQ t1 Þ þ kui kH1 ðXÞ Þ ¼ Oðui Þ; t0

t0

0

ð6:28Þ

sup kvðt; :ÞkH2 ðXÞ ¼ OðklvkH1 ðt0 ;t1 ;L2 ðXÞÞ þ kui kH2 ðXÞ Þ ¼ Oðui Þ

ð6:29Þ

sup kzc ðt; :ÞkH2 ðXÞ ¼ Oðkzc l þ cu2c kH1 ðt0 ;t1 ;L2 ðXÞÞ Þ; ¼ Oðui 2 Þ:

ð6:30Þ

t 0 6t6t 1

Proof. Let n0 6 1 be a smooth function in X, such that n0  1 in B2e and n0  0 in X n Be . Let n be the solution of

8 in Q; > < ot n ¼ DDn  Mn nðt; xÞ ¼ 0 on R; > : nð0; xÞ ¼ ðui =ui Þn0 ðxÞ in X:

Using (6.27) and (6.30), we get:

ð6:19Þ

Let us set r :¼ inf x2X1 nðT 0 ; xÞ. From the strong parabolic maximum principle, n > 0 in Q tt10 , and therefore, we get that r > 0, since X1 is a closed subset of X. Moreover, the parabolic maximum principle ~ in Q tt1 . In particular, we get v ~ðT 0 ; xÞ P ui r > 0 in also yields ui n 6 v 0 X1 . h qffiffiffiffi ðDþMÞeMt1 From Lemma 6.5, it follows that, for k P 2C , the term r s 2 2 2Mt1 1 2 0 2 e v ðx; T Þ  Cs k ðD þ MÞ ui e in (6.18) satisfies:

~ðx; T 0 Þ2  Cs1 k2 ðD þ MÞ2 ui 2 e2Mt1 P ui 2 v

r2 > 0 in X1 : 2

 Z t1 Z C 3=2 2 s k e2sg u3 y2 ui 2 t0 x  Z 0 þ e2sgðT ;xÞ ððDwðT 0 ;xÞÞ2 þ wðT 0 ;xÞ2 Þdx : ð6:20Þ

0

e2sgðT ;xÞ r2 dx 6

X1

Gl ðc; lÞ ¼ Oðui 4 Þ: Moreover, since u0 ¼ v, we have

~ Þ ¼ Gl ð0; l ~ Þ þ Gl ðc; lÞ þ 2hot zc ; ot v  ot v ~i 2 t 1 Gl ðc; l L ðQ Þ t0

0

~ðT 0 ; ÞiL2 ðXÞ þ 2hDzc ðT ; Þ; DvðT ; Þ  Dv

Using the fact that e2sg u3 remains bounded in obtain

krk2L2 ðX1 Þ 6

C ui 2

Z

t1

t0

Z

y2 þ

Z

x

we finally

ð6:21Þ ~, y ¼ ot ðv  v ~Þ, (6.21) implies the result of Recalling that w ¼ v  v Theorem 2.2. A.2. Proof of Theorem 2.4 Let uc be the solution of ðPl;c Þ, and let v be the solution of ðP l;0 Þ. Let us set zc :¼ v  uc . The function zc is a solution of

in Q; on R;

ð6:22Þ

in X:

It follows from the parabolic maximum principle that zc ðt; xÞ > 0 in Q tt10 . Thus uc 6 v in Q tt10 . Using the result of Lemma 6.4, we thus obtain:

uc ðt; xÞ 6 ui eMt1 in Q tt10 :

ð6:32Þ

~, and Thus, using (6.26) and (6.29), which are true for both v and v (6.27), (6.30), (6.31), together with Cauchy–Schwarz inequality, we get:

ð6:33Þ

A.3. Definitions of the state spaces E and of their neighbourhood systems In Examples 1 and 2 of Section 4, E is defined by

( E :¼

q 2 M; qðxÞ ¼

89 X

)

ak vk ðxÞ in X1 ; and qðxÞ ¼ m in X n X1 ;

k¼10

 ððDwðT 0 ; xÞÞ2 þ wðT 0 ; xÞ2 Þdx :

X

8 2 > < ot zc  DDzc ¼ zc lðxÞ þ cuc zc ðt; xÞ ¼ 0 > : zc ð0; xÞ ¼ 0

0

~ðT 0 ; ÞiL2 ðXÞ : þ 2hzc ðT 0 ; Þ; vðT 0 ; Þ  v

X

Q tt10 ,

ð6:31Þ

~ Þ  Gl ð0; l ~ Þj ¼ Oðui 3 Þ: jGl ðc; l

We deduce that, for k large enough

Z

t 0 6t6t 1

ð6:23Þ

k where vk are the characteristic functions of the intervals ð100 ; kþ1 Þ, 100 and ak are real numbers taken in finite subsets of ½m; M. In Example 1, ak 2 f1; 2g, and two distinct elements m1 , m2 of E, P P89 with m1 ¼ 89 k¼10 a1;k vk ðxÞ and m2 ¼ k¼10 a2;k vk ðxÞ in X1 , are defined as neighbours if and only if there exists a unique integer k0 in ½10; 89 such that a1;k0 6¼ a2;k0 . Note that, in this case, the number

of elements in E is 280 . In Example 2, ak 2 fðM  mÞj=20 þ m; for j ¼ 0 . . . 20:}, and the neighbourhood system is defined as follows: two distinct elements P P89 m1 , m2 of E, with m1 ¼ 89 k¼10 a1;k vk ðxÞ and m2 ¼ k¼10 a2;k vk ðxÞ in X1 , are neighbours if and only if (i) there exists a unique integer k0 in ½10; 89 such that a1;k0 6¼ a2;k0 , (ii) additionally, ja1;k0  a2;k0 j ¼ ðM  mÞ=20. Note that, in such a situation, the number of elements in E is 2180 . In Examples 3 and 4 of Section 4, E is defined by

( E :¼

q 2 M; qðxÞ ¼

17 X

)

ai;j vi;j ðxÞ in X1 ; and qðxÞ ¼ m in X n X1 ;

i;j¼2

Let k be the solution of

(

0

k ¼ kM þ cui 2 e2Mt1 on Rþ ;

ð6:24Þ

kð0Þ ¼ 0:

Since Kðt; xÞ :¼ kðtÞ is a supersolution of (6.22), we obtain that zc ðt; xÞ 6 kðtÞ in Q t01 . Thus, since k is increasing:

zc ðtÞ 6 kðt 1 Þ ¼

cui 2 e2Mt1 M

ðe

Mt1

 1Þ in

Q tt10 :

ð6:25Þ

Standard parabolic estimates (see e.g. [16]) then imply, using (6.15), (6.23), (6.25) and the hypothesis kui kC 2 ðXÞ 6 ui , that:

kot vkL2 ðQ t1 Þ ¼ OðklvkL2 ðQ t1 Þ þ kui kH1 ðXÞ Þ ¼ Oðui Þ; t0

0

t0

2

2

kot zc kL2 ðQ t1 Þ ¼ Oðkzc l þ cuc kL2 ðQ t1 Þ Þ ¼ Oðui Þ; t0

t0

ð6:26Þ ð6:27Þ

where vi;j are the characteristic functions of the square cells ði=20; i=20 þ 1=20Þ  ðj=20; j=20 þ 1=20Þ, and ai;j are real numbers taken in finite subsets of ½m; M. In Example 3, ai;j 2 f1; 2g. In this case, the number of elements in

m1

E is 2256 . Two distinct elements m1 , m2 of E, with P P17 2 1 ¼ 17 i;j¼2 ai;j vi;j ðxÞ and m2 ¼ i;j¼2 ai;j vi;j ðxÞ for x 2 X1 , are defined

as neighbours if and only if there exists a unique couple ði0 ; j0 Þ of integers comprised between 2 and 17 such that a1i0 ;j0 6¼ a2i0 ;j0 . In Example 4 ai;j 2 fðM  mÞj=20 þ m; for j ¼ 0 . . . 20:}. The number of elements in E is 21256 . In this case, two distinct elements P P17 2 1 m1 , m2 of E, with m1 ¼ 17 i;j¼2 ai;j vi;j ðxÞ and m2 ¼ i;j¼2 ai;j vi;j ðxÞ for x 2 X1 , are defined as neighbours if and only if (i) it exists a unique couple ði0 ; j0 Þ of integers comprised between 2 and 17 such that a1i0 ;j0 6¼ a2i0 ;j0 ; and (ii) additionally ja1i0 ;j0  a2i0 ;j0 j ¼ ðM  mÞ=20.

166

M. Cristofol, L. Roques / Mathematical Biosciences 215 (2008) 158–166

References [1] A.M. Liebhold, W.L. MacDonald, D. Bergdahl, V.C. Mastro, Invasion by exotic forest pests: a threat to forest ecosystems, For. Sci. Monogr. 30 (1995) 1. [2] National Invasive Species Information Center, USDA 2006, Executive Order 13112. [3] J.G. Skellam, Random dispersal in theoretical populations, Biometrika 38 (1951) 196. [4] N. Shigesada, K. Kawasaki, Biological Invasions: Theory and Practice, Oxford Series in Ecology and Evolution, Oxford University, Oxford, 1997. [5] P. Turchin, Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants, Sinauer Associates, Sunderland, MA, 1998. [6] A. Okubo, S.A. Levin, Diffusion and Ecological Problems – Modern Perspectives, second ed., Springer-Verlag, 2002. [7] R.A. Fisher, The wave of advance of advantageous genes, Ann. Eugenics 7 (1937) 355. [8] A.N. Kolmogorov, I.G. Petrovsky, N.S. Piskunov, Étude de l’équation de la diffusion avec croissance de la quantité de matière et son application à un problème biologique, Bull. Univ. État Moscou, Série Int. A 1 (1937) 1. [9] H. Berestycki, F. Hamel, L. Roques, Analysis of the periodically fragmented environment model. I. Species persistence, J. Math. Biol. 51 (1) (2005) 75. [10] H. Berestycki, F. Hamel, L. Roques, Analysis of the periodically fragmented environment model. II. Biological invasions and pulsating travelling fronts, J. Math. Pures Appl. 84 (8) (2005) 1101. [11] N. Kinezaki, K. Kawasaki, N. Shigesada, Spatial dynamics of invasion in sinusoidally varying environments, Popul. Ecol. 48 (2006) 263. [12] N. Shigesada, K. Kawasaki, E. Teramoto, Traveling periodic waves in heterogeneous environments, Theor. Popul. Biol. 30 (1986) 143.

[13] R.S. Cantrell, C. Cosner, Spatial ecology via reaction–diffusion equations, Series In Mathematical and Computational Biology, John Wiley and Sons, Chichester, Sussex, UK, 2003. [14] L. Roques, R.S. Stoica, Species persistence decreases with habitat fragmentation: an analysis in periodic stochastic environments, J. Math. Biol. 55 (2007) 189. [15] L. Roques, M.D. Chekroun, On population resilience to external perturbations, SIAM J. Appl. Math. 68 (1) (2007) 133. [16] L.C. Evans, Partial Differential Equations, University of California, Berkeley, AMS, 1998. [17] L. Roques, F. Hamel, Mathematical analysis of the optimal habitat configurations for species persistence, Math. Biosci. 210 (2007) 34. [18] B. Hajek, Cooling schedules for optimal annealing, Math. Oper. Res. 13 (1988) 311. [19] D. Henderson, S.H. Jacobson, A.W. Johnson, The Theory and Practice of Simulated Annealing, in: F. Glover, G. Kochenberger (Eds.), Handbook on Metaheuristics, Kluwer Academic Publishers, Norwell, MA, 2003, p. 287. [20] S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671. [21] P. Albano, D. Tataru, Carleman estimates and boundary observability for a coupled parabolic–hyperbolic system, Electron. J. Different. Eqn. 22 (2000) 1. [22] A. Fursikov, Optimal control of distributed systems, Translations of Mathematical Monographs, American Mathematical Society, Providence, RI, 2000. [23] L. Baudouin, J.P. Puel, Uniqueness and stability in an inverse problem for the Schrödinger equation, Inv. Prob. 18 (2002) 1537. [24] M. Cristofol, P. Gaitan, H. Ramoul, Inverse problems for a two by two reaction– diffusion system using a Carleman estimate with one observation, Inv. Prob. 22 (2006) 1561. [25] L. Cardoulis, M. Cristofol, P. Gaitan, Inverse problem for a Schrödinger operator in an unbounded strip, J. Inv. Ill-Posed Prob. V (16) (2008) 127.