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

1 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.

2 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₁ and x₂ is ∫x1x2u(x,t)dx . Therefore, the total number of individuals at any time t is ∫m0Mu(x,t)dx , where m₀ 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₁ to t₂ the number ∫t1t2∫x1x2m(x)u(x,t)dx of individuals having mass between x₁ and x₂ 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 = 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₁ to t₂ the total number of migrants with masses between x₁ and x₂ is ∫t1t2∫x1x2f(x,t)dx .

Thus, we obtain the conservation law

for all (x₁, t₁), (x₂, t₂) ∈ Ω ≔ {(x,t) : m₀ ≤ x ≤ M, 0 ≤ t < +∞).

Assuming smoothness of u and g, as well as continuity of m and f, equation (1) may be transformed into the single PDE

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₁ and time t₂ is ∫t1t2∫m0Mb(x,t)u(x,t)dxdt . Since this quantity must also be ∫t1t2g(m0)u(m0,t)dt , we obtain the renewal condition

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]).

3 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

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 ↦ g(x)⁻¹ is Lebesgue integrable in (m₀, M). Therefore, the time required for individuals to grow from mass x₁ to mass x₂ is ∫x1x2g(x)−1dx . Characteristics are the fundamental concept in the analysis of hyperbolic problems because PDEs simplify to ODEs along these curves. Denote G(x)=∫m0xg(s)−1ds . Since the function G is strictly increasing, and thus it is invertible, the characteristic curve, in coordinate system, passing through the point (x,t) ∈ Ω is easily found to be

Here and subsequently, for notational simplicity, we drop the variables x, t and write ξ(τ) instead of ξ(τ; x,t) but keeping in mind that ξ(τ) depends on the choice of (x,t).

Differentiating the solution u along the characteristics yields

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

where u(ξ(τ₀), r₀) is to be determined.

Make the change of variables ξ₁ = ξ(τ₁), g(ξ₁)⁻¹ dξ₁ = dτ₁ to obtain

and therefore

Denote π(x):=exp−∫m0xm(ξ)g(ξ)−1dξ , which is the probability to grow from mass mo, which is the mass of all newborns, to mass x without taking into account harvesting. Then π(x)π(x0)−1=exp−∫x0xm(ξ)g(ξ)−1dξ , for any x₀ < x, is the probability that an individual of mass x₀ will grow to mass x. We now assume that ∫m0Mm(x)g(x)−1dx=+∞ , while the function x ↦ m(x)g(x)⁻¹ is locally integrable in [m₀, M); that is, the probability to accumulate the maximum possible mass M equals zero. Thus π(M) = 0, and in the sequel, for simplicity of our investigation, we assume that π(x) = 0, and therefore u(x, t) = 0 for all x ≥ M, t > 0. Also, the product of n 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, τ = t, that is, considering the previous equation on the characteristics that emanate from points (ξ(0),0) on the x axis, where ξ(0) = G⁻¹ (G(x) – t), gives

Similarly, taking τ₀ = t – G(x), τ = t, that is, considering equation (9) on the characteristics that emanate from points (m₀, t – G(X)), gives

4 Reduction of the inverse problem to a system of integral equations

Denote the number of births in unit time by B(t) (i.e., B(t) ≔ g(m₀)u(m₀, t)). Then substituting B(t) and condition (4) into equations (10), (11) yields

Now, substitute the last two equations into condition (3) to obtain

Similarly, substituting equations (12), (13) into condition (6) gives

Note that the first term ∫m0G−1(t)w(x)u(x,t)dx represents the number (if w(x) ≡ 1) or the biomass (if w(x) ≡ x)of individuals born between time zero and time t and surviving to time t, and the second term ∫G−1(t)+∞w(x)u(x,t)dx represents the number of individuals born before the time t = 0 and still surviving at time t. Since the time required for individuals to grow from birth to the maximum possible mass M is G(M)=∫m0Mg(x)−1dx , and the probability to accumulate the maximum possible mass M equals zero (because G(M) is the maximum lifetime), we see that the integral ∫G−1(t)+∞w(x)u(x,t)dx equals zero whenever t ≥ G(M). This fact also follows from our convention that π(x) = 0, and therefore u(x, t) = 0, for all x ≥ M, t > 0.

Note also that exp (−∫0tE(τ1)dτ1) ) is the probability that an individual will not be harvested up to time t. Using the notation E(t):=exp⁡∫0tE(τ1)dτ1 and multiplying equation (14) by 𝓔 (t), we rewrite this equation in the form

Denoting h(τ2;x,t):=f(ξ(τ2),τ2)g(ξ(τ2))g(x)π(x)π(ξ(τ2)), q1(t):=∫G−1(t)+∞b(x,t)u0(ξ(0))g(ξ(0))g(x)π(x)π(ξ(0))dx , and also 𝓑(t) ≔ B(t)𝓔(t), we simplify the last equation and obtain

Changing the order of integration and making the change of variables τ = t – G(x), dτ = –g(x)⁻¹ dx in the first integral give

where we take k₁ (t, τ) ≔ b(G⁻¹(t – τ), t)π(G⁻¹(t – τ)), k2(t,τ):=∫G−1(t−τ)+∞b(x,t)h(τ;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(τ; x, t) = 0 for all x ≥ M, τ < t, thus we may replace all the improper integrals by proper ones.

Similarly we deal with equation (15) to obtain

where we use the notation k^3(t,τ):=w(G−1(t−τ))π(G−1(t−τ)), k^4(t,τ):=∫G−1(t−τ)+∞w(x)h(τ;x,t)dx, and q^2(t):=∫G−1(t)+∞w(x)u0(ξ(0))g(ξ(0))g(x)π(x)π(ξ(0))dx . Dividing the last equation by p(t) gives

where k3(t,τ):=k^3(t,τ)p(t)−1,k4(t,τ):=k^4(t,τ)p(t)−1,q2(t):=q^2(t)p(t)−1. .

5 Existence-uniqueness theorem

We formulate our main result as a theorem.