Inverse problem for a physiologically structured population model with variable-effort harvesting

Abstract We consider the inverse problem of determining how the physiological structure of a harvested population evolves in time, and of finding the time-dependent effort to be expended in harvesting, so that the weighted integral of the density, which may be, for example, the total number of individuals or the total biomass, has prescribed dynamics. We give conditions for the existence of a unique, global, weak solution to the problem. Our investigation is carried out using the method of characteristics and a generalization of the Banach fixed-point theorem.


Introduction
Population dynamics has traditionally been the dominant branch of mathematical biology. Population models play a critical role in helping us to understand the dynamic processes involved, in making practical predictions, and thus in better understanding the natural world. The mathematics we use to study such models mainly involves evolution equations. Therefore, ODE and PDE problems describing the behavior of biological systems that depend on a continuous time variable are of great importance in the natural sciences, and they have been extensively studied using a variety of mathematical techniques. The benefit of PDE models is that they can take into account the structure of populations, provided that this structure influences population size and growth in a major way.
The central problem in population ecology is to determine how the age structure of a population evolves in time and to understand the factors that regulate animal and plant populations. But age is just one of the many structure variables that demographers and ecologists study. Any quantity that is characteristic of an individual's state can be used as a structural variable; these include maturation level, size, weight, and so on. Thus age-structured models are akin to more general physiologically structured models. For example, in many non-mammalian populations the evolution of the population, especially the mortality rate, certainly depends upon the size of the animals; the survival probability is small for young fish or insects. Of course, the study of age-structured models is considerably simpler than the study of general size-structured or mass-structured models, primarily because age increases linearly with the passage of time while the linkage of size or mass with time may be less predictable. The use of physiologically structured models to describe biological systems has attracted the interest of many researchers and has a long standing tradition. The books by Charlesworth [1], Metz and Diekmann [2], Diekmann [3], Cushing [4], Kot [5], and Murray [6] give a good survey as well as the wide spectrum of applicability of such models.
In many real-life problems all of the input (functional) parameters may not be known a priori and cannot be directly observed. So in parameter identification problems (a subclass of inverse problems) we ask whether it is possible to take certain additional measurements or, mathematically, impose an additional condition that describes the modeling process and thereby determine unknown parameters from additional experimental data. Inverse problems are currently of great research interest, and thus the development of efficient analytical and numerical methods to deal with inverse problems in population biology is an important task of applied mathematics. The direct problems of structured populations have been extensively studied for many kinds of models using a variety of mathematical techniques. But in recent years many researchers have focused their attention on developing methodologies for solving inverse problems governed by structured population models (see, e.g., Ackleh [7], Banks [8], Pilant and Rundell [9], Rundell [10], Gyllenberg et al. [11], Perthame and Zubelli [12], Wood and Nisbet [13]).
Thus our research is devoted to the inverse problem of determining how the physiological structure of a harvested population evolves in time, and of finding the time-dependent effort to be expended in harvesting, so that the weighted integral of the density, which may be, for example, the total number of individuals or the total biomass, has prescribed dynamics. Note that we study the simple case when the model equation is a single first-order PDE; such an equation is hyperbolic. And the fundamental idea associated with hyperbolic equations is the notion of a characteristic, a curve in space-time along which signals propagate. So, using the method of characteristics, which is highly effective for investigating hyperbolic continuous-time models (see, e.g., Logan [14], Brauer and Castillo-Chavez [15]), the solution of the inverse problem can be expressed as a fixed point of some appropriately chosen integral operator in a suitable metric space. Then the Banach fixed-point theorem being an important tool in the theory of metric spaces is used to show the existence of a unique solution to the original inverse problem. Furthermore, this methodology allows us to generate a sequence of iterates, called the Picard iterates, that under certain conditions converges to the solution of the inverse problem.

Population model with physiological distribution
Suppose that a population is structured by mass (or some other physiological quantity, say, maturation level, size, and so on). Any quantity that is characteristic of an individual's state can be used as a structural variable. Let u.x; t/ be the unknown density of the population at time t with respect to a mass variable x, so that the number of individuals at time t having mass between x 1 and x 2 is R x 2 x 1 u.x; t / dx. Therefore, the total number of individuals at any time t is R M m 0 u.x; t / dx, where m 0 is the mass of all newborns and M is the maximum possible mass that individuals can accumulate over their lifespan. Assume that mass is accumulated by individuals of mass x with the growth rate g.x/ (i.e., all individuals of the same mass experience the same growth rate). We also assume that members leave the population through death, and that there is an mass-dependent per capita death rate m.x/. This means that over the time interval from t 1 to t 2 the number R t 2 t 1 R x 2 x 1 m.x/u.x; t / dx of individuals having mass between x 1 and x 2 die. But a population may be changed by migration either into or out of the population. We assume that the rate of migration is described by the function f D f .x; t /, with a positive value of f denoting immigration and a negative one denoting emigration; that is, over the time interval from t 1 to t 2 the total number of migrants with masses between x 1 and x 2 is

Thus, we obtain the conservation law
for all .x 1 ; t 1 /; .
Assuming smoothness of u and g, as well as continuity of m and f , equation (1) may be transformed into the single PDE @ @t u.x; t / C @ @x g.x/u.x; t / Á D m.x/u.x; t / C f .x; t /; .x; t / 2 : Next, we assume that the birth process is governed by a function b.x; t / called the per capita birth rate. Thus, the total number of births between time t 1 and time t 2 is R t 2 t 1 R M m 0 b.x; t /u.x; t / dxdt. Since this quantity must also be R t 2 t 1 g.m 0 /u.m 0 ; t/ dt , we obtain the renewal condition b.x; t /u.x; t / dx; t 0: In order to complete the model, we must specify an initial age distribution In summary, the PDE (2) and the two auxiliary conditions (3), (4) define a well-posed mathematical problem to determine how the mass structure of the entire population evolves (see [6,14,15]).

Formulation of an inverse problem. Definition of a weak solution
We wish to study the effect on the model presented above of the removal of members of the population. Now suppose that the model is subject to a time-dependent harvesting at the per capita rate E.t /. Then the harvested population is modeled by the PDE along with conditions (3), (4). This type of harvesting arises in the modeling of fisheries, where it is often assumed that the number of fish caught per unit time is proportional to E.t /, the effort expended in fishing. This fishing effort may be measured, for example, by the number of boats fishing at a given time. The study of harvesting in population models can be found in [6,15]. Now we consider the inverse problem of determining how the physiological (e.g., mass) structure of a harvested population evolves in time, and of finding the time-dependent effort to be expended in harvesting, so that the weighted integral of the density, which may represent the total number of individuals or the total biomass, has prescribed dynamics given by Z M m 0 w.x/u.x; t / dx D p.t /; t 0: That is, in the inverse problem (the parameter identification problem) we seek both the mass density u and the harvest rate E that satisfy the PDE (5), along with conditions (3), (4), (6). Note that the physiological variable and time are related by the characteristic equation This equation has the set of solution curves (called characteristics) where g is assumed to be a positive function such that x 7 ! g.x/ 1 is Lebesgue integrable in .m 0 ; M /. Therefore, the time required for individuals to grow from mass x 1 to mass x 2 is R x 2 x 1 g.x/ 1 dx. Characteristics are the fundamental concept in the analysis of hyperbolic problems because PDEs simplify to ODEs along these curves. Denote G.x/ D R x m 0 g.s/ 1 ds. Since the function G is strictly increasing, and thus it is invertible, the characteristic curve, in coordinate system, passing through the point .x; t / 2 is easily found to be Here and subsequently, for notational simplicity, we drop the variables x, t and write . / instead of . I x; t / but keeping in mind that . / depends on the choice of .x; t /.
Differentiating the solution u along the characteristics yields d d u. . /; / D @ @t u. . /; / C g. . // @ @x u. . /; /: Then, using the previous relation, equation (5) can be transformed to Then, using the variation of parameters method, which works for all linear equations, we can find a particular solution of (8) and thus the general solution to (8) has the form and therefore , which is the probability to grow from mass m 0 , which is the mass of all newborns, to mass x without taking into account harvesting. Then .
for any x 0 < x, is the probability that an individual of mass x 0 will grow to mass x. We now assume that R M m 0 m.x/g.x/ 1 dx D C1, while the function x 7 ! m.x/g.x/ 1 is locally integrable in OEm 0 ; M /; that is, the probability to accumulate the maximum possible mass M equals zero. Thus .M / D 0, and in the sequel, for simplicity of our investigation, we assume that .x/ D 0, and therefore u.x; t / D 0 for all x M , t > 0. Also, the product of and any other function is interpreted to be zero whenever x M (even if the latter function is not defined on this interval).
Taking 0 D 0, D t, that is, considering the previous equation on the characteristics that emanate from points x m 0 ; 0 Ä t Ä G.x/: (10) Similarly, taking 0 D t G.x/, D t , that is, considering equation (9) on the characteristics that emanate from points .m 0 ; t G.x//, gives x m 0 ; G.x/ Ä t < C1: (11) x m 0 ; G.x/ Ä t < C1: (13) Now, substitute the last two equations into condition (3) to obtain Similarly, substituting equations (12), (13) into condition (6) gives t 0: (15) Note that the first term Note also that exp is the probability that an individual will not be harvested up to time t.
Using the notation E.t / WD exp R t 0 E. 1 / d 1 and multiplying equation (14) by E.t /, we rewrite this equation in the form Changing the order of integration and making the change of variables D t G.x/, d D g.x/ 1 dx in the first integral give where we take x; t /h. I x; t / dx. We recall that the product of and any other function is interpreted to be zero whenever x M (even if the latter function is not defined on this interval). Therefore, for example, h. I x; t / D 0 for all x M , < t, thus we may replace all the improper integrals by proper ones.
Similarly we deal with equation (15) (12), (13). Note that the function E is required to be at least positive locally absolutely continuous in OE0; C1/ (suitable assumptions on given functions guarantee this), then E is differentiable a.e. in this interval and E 0 is locally Lebesgue integrable (see [16]).

Existence-uniqueness theorem
We formulate our main result as a theorem. We remark that any function with continuous first derivatives is locally Lipschitz. Thus the property of being locally Lipschitz is stronger than continuity, yet weaker than continuous differentiability. where k k 1;T is the uniform norm, that is, and k k L;T is the Lipschitz constant, which serves as a seminorm, that is, Let us show that the mapping P has precisely one fixed point in the Banach space introduced above. To this end, by a generalization of the Banach fixed-point theorem, we must show that the considered space is invariant under the mapping P, and some its iterate P n is a contraction. It is easily seen that the functions P 1 .B; E/ and P 2 .B; E/ are nonnegative and Lipschitz in OE0; T if so are the functions B and E; that is, P.B; E/ 2 M.T / whenever .B; E/ 2 M.T /. Further, for any pairs .B 1 ; E 1 /, .B 2 ; E 2 / in M.T /, we obtain the estimates and kP i .
for all t 2 OE0; T and every i 2 f1; 2g. Here we use the notation ƒ. : Therefore, we have kP.
Á d for all t 2 OE0; T , where constants C 1 , C 2 are easily determined. From the last estimate it follows that for all t 2 OE0; T . Thus there exists a positive constant C 3 not depending on on the choice of .B 1 ; E 1 / and .B 2 ; E 2 / such that for all t 2 OE0; T . From this estimate we obtain and therefore and so on. Mathematical induction can be used to prove that the following estimate holds for all positive integers n and all t 2 OE0; T : Thus, in general, the mapping P itself is not a contraction on the considered Banach space, but some its iterate P 2n is a contraction provided that n is chosen sufficiently large to satisfy the inequality .C 3 T / n nŠ < 1. Consequently, for all positive T , using a generalization of the Banach fixed-point theorem, the mapping P has precisely one fixed point in the Banach space M.T /. This fixed point is a solution of the system (16), (17) in the interval OE0; T . Taking into account that the value T may be chosen arbitrarily large, we may extend this solution to any finite time T 1 , where T 1 > T , and therefore we obtain a global existence result, guaranteeing the existence of a unique solution .B; E/ of the system (16), (17) for all t in OE0; C1/, where B and E are nonnegative locally Lipschitz. Then applying proposition 4.1 gives the unique, global, weak solution .u; E/ to the inverse problem (3)-(6), where u is nonnegative locally Lipschitz and E is locally essentially bounded in their respective domains. The latter follows from the Lipschitz property of E in each interval OE0; T for all T > 0, and therefore the fact that E 0 is essentially bounded in each such interval. In fact, for any open bounded interval I , a function belongs to the Sobolev space W 1;1 .I / if and only if it admits a Lipschitz continuous representative. In particular, this representative is differentiable a.e. in I and its derivative is essentially bounded as an element of L 1 .I / (see [16]). Thus the theorem is proved.
Next we establish the validity of the conservation law for the weak solution to the inverse problem.
Proposition 5.2. The weak solution whose existence is guaranteed by theorem 5.1 also satisfies the integral form of the conservation law f .x; t / dx dt (18) for all .x 1 ; t 1 /; .x 2 ; t 2 / 2 .
We now integrate equation (5) over f .x; t / dx dt: Note that the possibility of integrating follows from the essential boundedness of both sides of (5)  It is clear that some additional restriction must be placed on the function p in order to guarantee the positivity of E.t/ for all t . Since we proved the existence of a solution to the system (16), (17) without exhibiting it explicitly, this restriction on p cannot be determined in an explicit form as well. But in simple cases when the solution .B; E/ of the system (16), (17) can be determined explicitly, we can write a sufficient and necessary condition on p that the function E is nonnegative or, equivalently, E 0 .t /E.t / 1 0 for all t . Let us consider, for example, the simple case when the functions b, m are constant, that is, b.x; t / Á b, m.x/ Á m (in this case we formally put M D C1). Also assume, as in an age-structured model, that g.x/ Á 1, w.x/ Á 1, m 0 D 0, and f .x; t / Á 0; the latter assumption means that there is no migration. Then equation (16) can be rewritten as and hence be reduced to the ordinary differential equation Remark 5.4 (Approximation of the solution). The Banach fixed-point theorem, as well as its generalization used above, provide an approximation method for solving the system of integral equations (16), (17) and therefore the inverse problem (3)-(6).
Precisely, for any nonnegative Lipschitz functions B 0 and E 0 , which is taken to be an initial approximation, the sequence of iterates .B 0 ; E 0 /; P.B 0 ; E 0 /; P.P.B 0 ; E 0 //; : : : converges in the appropriate norm topology to the fixed point of P. Thus we may define the iteration scheme for all t 0 and nonnegative integers n: Proceeding in this manner, we generate the sequence .B 1 ; E 1 /; .B 2 ; E 2 /; .B 3 ; E 3 /; : : : of iterates, called the Picard iterates, that, under conditions of theorem 5.1, converges to the solution of the system (16), (17) in the k k M.T / norm for all T > 0. Denoting this solution by .B ; E /, we have k.B n ; E n / .B ; E /k M.T / ! 0 as n ! C1, and therefore kB n B k 1;T , kB n B k L;T , as well as kE n E k 1;T , kE n E k L;T , approach zero as n ! C1. The sequence .B n ; E n / obtained above gives the sequence of approximate solutions .u n ; E n / to the inverse problem (3)-(6); the sequence u n is determined explicitly by the right sides of (12), (13), where the functions B and E must be replaced by B n and E n , respectively; B n is found by the formula B n .t / D B n .t /E n .t / 1 , and E n is taken to be E n .t / D E 0 n .t /E n .t / 1 (see proposition 4.1). Since E n , for all positive integers n, is Lipschitz continuous in OE0; T , which follows from its representation, we see that E 0 n belongs to L 1 .0; T /; that is, E 0 n is essentially bounded in .0; T /. Let .u ; E / denote the exact weak solution to the inverse problem. By the convergence of iterates, we conclude that kE 0 n E 0 k L 1 .0;T / , kB n B k 1;T approach zero as n ! C1, and consequently, kE n E k L 1 .0;T / , ku n u k 1; .T / both go to zero for all T > 0.
To summarize, we have constructed the sequence of approximate solutions .u n ; E n / to the inverse problem (3)- (6), and this sequence tends to the exact weak solution .u ; E / in the appropriate norm topology.
Note that in equation (5) the death term is assumed to be linear. Such an assumption appears to be reasonable for simple organisms such as micro-organisms. For more complicated organisms like plants, animals, or human beings this is obviously an oversimplification since it ignores intra-species competition for resources. (25) where the logistic term n R M m 0 w.s/u.s; t / ds Á provides a mechanism for increased mortality as the population increases, which may be explained by intra-species competition for resources.
The last equation combined with condition (6) gives the linear equation of the form @ @t u.x; t / C @ @x g.x/u.