Stability and Optimization in Structured Population Models on Graphs

We prove existence and uniqueness of solutions, continuous dependence from the initial datum and stability with respect to the boundary condition in a class of initial--boundary value problems for systems of balance laws. The particular choice of the boundary condition allows to comprehend models with very different structures. In particular, we consider a juvenile-adult model, the problem of the optimal mating ratio and a model for the optimal management of biological resources. The stability result obtained allows to tackle various optimal management/control problems, providing sufficient conditions for the existence of optimal choices/controls.


Introduction
This paper is devoted to the following initial-boundary value problem for a system of balance laws in one space dimension: x) ∈ R + × R + g i (t, 0) u i (t, 0+) = B i t, u 1 (t), . . . , u n (t) t ∈ R + u i (0, x) = u o i (x) x ∈ R + i = 1, . . . , n . (1.1) Here, i = 1, . . . , n and t ∈ R + is time. The "space" variable x varies in R + and in the applications of (1.1) will have the meaning of a biological age, or size. The scalar functions g 1 , . . . , g n are growth functions, d 1 , . . . , d n are the death rates and u o 1 , . . . , u o n constitute the initial data. A key role is played by our choice of the birth function B i , for i = 1, . . . , n, which we assume of the form B i (t, u 1 , . . . , u n ) = α i t, u 1 (x 1 −), . . . , u n (x n −) + β i I1 u 1 (x) dx , . . . , In u n (x) dx (1.2) for suitable functions α i , β i , pointsx i > 0 and measurable I i ⊆ R + , for i = 1, . . . , n.
The literature on equations similar to (1.1) is vast. We refer for instance to the exhaustive monograph [14] or to the more recent edition of [5] and to the references therein. Specific features of (1.1) are that it is a system, boundary conditions may contain both a local term, the α i , and a nonlocal term, the β i .
From the analytical point of view, in the present treatment we emphasize the role of the total variation, setting the main result in BV. In particular, this allows to consider a function of the type (1.2) and to prove that the boundary data are attained in the sense of traces, also due to the boundary being non characteristic. In this setting, the stability of solutions with respect to α i and β i is also obtained. Moreover, the techniques used in the sequel can easily be extended to more general source terms as well as to situations where also the space distribution needs to be taken into account.
From the modeling point of view, the use of boundary conditions of the type (1.2) unifies the treatment of rather diverse situations. First, it comprises the standard case always covered in the literature on renewal equations, where the independent variable x varies along a segment or a half line, see Figure 1, left. The dependent variable u represents the population density that at time t is of size (or age) x. A more complicate structure was recently considered in [1], see Figure 1, right. There, the size/age biological variable varies along a graph consisting of 2 distinct sets, corresponding to the juvenile and to the adult stages in the development of the considered species. Here, we are able to deal also with this situation, as depicted in Figure 2, right. . Left, a framework corresponding, for instance, to sexual reproduction: the two branches correspond to males M and females F . Right, a structure possibly accounting for the exploitation of biological resources: when juveniles reach the adult stage, they are split into a part S which is, say, sold and a part R used for reproduction.
The finite propagation speed intrinsic to models of the type (1.1) clearly allows to combine various instances of the graphs above. Other situations of biological interest can be for instance a three stage linear structure or a tree shaped one, see Figure 3. These schemes, as well as many others, all fit into the scope of Theorem 2.4 below. In this connection, we recall that similar network structures are widely considered in the framework of vehicular traffic modeling, see [9].
In the case of nonlinear systems of balance laws with flow independent from the space variable, the initial boundary value problem has been widely investigated, see for instance [7]. For the relations between the problems with boundaries and with junction see [8,Proposition 4.2].
The present treatment is self-contained. Section 2 is devoted to the analytical results. Specific applications are in Section 3, where sample numerical integrations are also provided. All technical details are deferred to Section 4.

Analytical Results
Throughout, we use the standard notation R + = [0, +∞[ andR + = ]0, +∞[. When A and B are suitable subsets of R m , C 0 (A; B), respectively C 0,1 (A; B), L 1 (A; B) or L ∞ (A; B), is the set of continuous, respectively Lipschitz continuous, Lebesgue integrable or essentially bounded, maps defined on A and attaining values in B. For the basic theory of BV functions we refer to [3].
When referring to a function u : R + × R + → R, the first argument is time, the second is the biological age/size variable. If I ⊆ R + is an interval, we denote Preliminarily, we consider the following initial-boundary value problem for a linear scalar balance law, or renewal equation in [14,Chapter 3]: under the following assumptions for positiveǧ,ĝ and sup t∈R + TV g(t, ·) < +∞ sup t∈R + TV ∂ x g(t, ·) < +∞ ; The solutions to (2.1) can be written in terms of the ordinary differential equationẋ = g(t, x). If g satisfies (g), we can introduce the globally defined maps Denote γ(t) = X(t; 0, 0), its inverse being t = Γ(x). Note that Recall the following definition of solution to (2.1), see also [4,6,11,14,18].
The following Lemma summarizes various properties of the solution to (2.1), see also [14]. Here, we stress the role of BV estimates. The proof is deferred to Section 4.
solves (2.1) in the sense of Definition 2.1. Moreover, there exists a constant C dependent only on g and d, see (4.6), such that the following a priori estimates hold for all t ∈ R + : x)];R) e Ct x < γ(t) . (2.7) Moreover, for any interval I ⊆ R + , For every t ∈ R + , there exists a positive L dependent onǧ,ĝ, C and For u o , u o ∈ (L 1 ∩ BV)(R + ; R) and b , b as in (b), the solutions u and u to satisfy the stability and monotonicity estimates Recall that in the present case of a linear conservation law, the definition of weak solution at 2. is equivalent to the definition of Kružkov solution [11,Definition 1].
It is immediate to verify that for u o = 0 and b = 0, problem (2.3) admits the solution u = 0. Hence, the monotonicity property (2.12) also ensures that non-negative initial and boundary data in (1.1)-(1.2) lead to non-negative solutions.
In order to pass to system (1.1), we need the following notation for norms and total variations of functions attaining values in R n : As a reference for the usual definition of weak solution to scalar conservation laws, see [4,11].
. . , g n satisfying assumptions (g) and in the sense of Definition 2.1.
The following result ensures the well posedness of (1.1)-(1.2). Its proof is presented in Section 4.
We now state separately the stability of solutions to (1.1)-(1.2) with respect to the birth function B. This result plays a key role in the optimization problems considered below.
satisfy the assumptions of Theorem 2.4. Then, the corresponding solutions u and u are such that The proof is deferred to Section 4. In applications of Theorem 2.4 to systems motivated by, for instance, structured population biology, further assumptions are natural and lead to further reasonable properties.
The proof follows immediately from Theorem 2.4 and from (2.12), hence it is omitted.

Applications
This section is devoted to sample applications of Theorem 2.4 and Theorem 2.5 to models inspired by structured population biology. We selected three cases corresponding to three different graphs, namely those in Figure 1, right, and in Figure 2. First, the well posedness ensured by Theorem 2.4 provides a ground for the reliability of each model. Then, the stability result in Theorem 2.5 allows to consider further problems. On the one hand, it ensures the existence of a choice of parameters in the equations that lead to solution that best approximate a given set of data. On the other hand, it allows to tackle the problem of optimal mating ratio in a population with sexual reproduction. Finally, we consider the problem of the optimal management of a biological resource. In the former case, the presentation is based on [1,2] where a sensitivity analysis for a model belonging to the class (1.1)-(1.2) is proved. In the latter cases, we provide numerical integrations showing further qualitative properties of the models considered.

A Nonautonomous Juvenile-Adult Model
In (1.1)-(1.2), setting n = 2 and with reference to the structure in Figure 1, right, in the case β = 1, which we state here for completeness: Theorem 2.4 then applies and ensures the well posedness of (3.2) under assumptions slightly different from those in [1].
Then, problem (3.2) admits a unique solution in the sense of Definition 2.3, the continuous dependence estimates (2.15)-(2.16) and the stability estimate (2.18) apply.
For completeness, we remark that the model in [1] contains the following slightly more general boundary inflow: still allows to apply Theorem 2.4. Indeed, with this variable, the second equation in (3.2) becomes which is again of the type (2.1) and hence Theorem 2.4 can still be applied. The stability proved above allows to tackle the problem of parameter identification. Indeed, through a Weierstraß argument based on Theorem 2.5, one can prove the existence of a set of parameters in (3.2) that minimizes a continuous functional representing the distance between the computed solution and a set of experimental data. For a detailed sensitivity analysis for a juvenile-adult model we refer to [2].

Optimal Mating Ratio
Consider a species consisting of males and females, whose densities at time t and age a are described through the functions M = M (t, a) and F = F (t, a) on a structure as that in Figure 2, left. A natural model is then Here, κ µ, respectively (1 − κ) µ, is the mortality rate of males, respectively females, with µ > 0 and κ ∈ [0, 1]. The positive parameter η ∈ [0, 1] defines the ratio of male to female newborns, in the sense that every η M males, (1−η) F females are born. The constant ν is the fertility rate. We describe the mating ratio at age a through the parameter ϑ, with ϑ ∈ [0, 1] as follows. The fertile ages are those in the intervals [m 1 , m 2 ] for males and [f 1 , f 2 ] for females, where m 1 , m 2 , f 1 , f 2 are positive constants. According to (3.4), all individuals in their fertile age might contribute to reproduction provided the condition imposed by the presence of the mating ratio ϑ is met. The proof is immediate and, hence, omitted. Here, we note that the presence of C 1 positive weights in the integrands defining the boundary data can be recovered through a change of variables entirely similar to that in (3.3).
A first immediate property of the solutions to (3.4) is that a zero initial density in either of the two sexes leads to the extinction of the other at exponential speed.
Several different optimization problems can be tackled in the framework of (3.4). It is possible to investigate the relations between the parameters κ (identifying relative mortality), η (the relative natality) and ϑ (the mating ratio). Below, we look for the optimal mating ratio for given relative natality and mortality coefficients.
To this aim, consider the instantaneous average fertility rate over the fertile population Remarkably, this leads to the maximal fertility rate coherently with the classical harmonic mean law, see [10,13,15,16,17,19].
On the other hand, the right hand sides in (3.5) and (3.6) are time dependent and it can be hardly accepted that ϑ is instantaneously adjusted to the value that maximizes R. More reasonably, one may imagine that ϑ is optimal 1 over a suitably long time interval. We are thus lead to introduce the utility function We thus consider the problem find ϑ that maximizes R(ϑ; T, M o , F o ) .
A straightforward corollary of Theorem 2.5 ensures the existence of one such ϑ.

Corollary 3.3.
Under the assumptions of Corollary 3.2, for any T ∈R + and any initial datum The graph of the average fertility rate R(ϑ; T, M o , F o ) as a function of ϑ for T = 500 is in Figure 4. The outcome shows a reasonable qualitative behavior. As ϑ → 0 or ϑ → 1, the number of newborns goes to 0; hence the population extinguish. Near to ϑ = 0.77 there is an optimal choice for the parameter ϑ with respect to the average fertility rate (3.7), which yields a maximal value of 0.83, see Figure 5 For completeness, we precise that the numerical integration above was obtained using a Lax-Friedrichs algorithm, see [12, § 12.5], with space mesh ∆a = 0.04167.

Management of a Biological Resource
In biological resource management, one typically rears/breeds a species up to a suitable stage, then part of the population is sold and part is used for reproduction. The equations (1.1)-(1.2) comprehend this situation. Indeed, call J = J(t, a) the density of the juveniles at time t of age or size a. Juveniles reaching the age/sizeā are then selected. The density S = S(t, a) refers to those individuals that are going to be sold, while R = R(t, a) stands for the density of those reserved for reproduction purposes. One is thus lead to the following model, defined on the structure in Figure 2, right:  Above, we used the obvious notation for the growth and mortality functions g J , g S , g R and d J , d S , d R . The birth rate is described through the function β. A key role is played by the parameter η ∈ [0, 1] which quantifies the percentage of juveniles selected for the market. System (3.10) fits into (1.1)-(1.2) setting Let β ∈ C 0,1 (R + ; R) be such that β(0) = 0. For any η ∈ [0, 1] and any initial data  Here, C J (a) is the unit cost to grow a juvenile at age a, and similarly C R , C S are the costs for the other two groups. The gain obtained selling an adult at age a is G(a). We denoted by J = J(t, a), S = S(t, a) and R = R(t, a) the solution to (3.10) with initial datum J o and S o = 0, R o = 0 and with the selection parameter η. A direct consequence of Theorem 2.5 is the following corollary. The proof relies on Weierstraß Theorem, exactly as that of Corollary 3.3 and is here omitted.
As a specific example, we consider the situation identified by the following choices of functions and parameters in (3.10)-(3.11): The graph of the cost G(η; T ) − C(η; T ) (see (3.11)) for T = 15 with respect to λ is in Figure 6. The outcome shows a reasonable qualitative behavior. As η → 0, nothing is sold, all population members are kept for reproduction, the population increases exponentially as also does the functional G(η; T ) − C(η; T ). On the contrary, for η → 1, all population members are immediately sold giving a positive gain and the population vanishes as also G − C. Near to η ≈ 0.23 there is an optimal balance, given the chosen unitary costs and gain (3.11)-(3.12).
With the chosen parameters, the optimal choice for η is η * ≈ 0.23, which yields a gain of about 260.48 at time t = 15. As expected, different choices of η have deep influences on the solutions to (3.10)-(3.12)-(3.13), as shown in Figure 7.
For completeness, we precise that the numerical integration above was obtained using a Lax-Friedrichs algorithm, see [12, § 12.5], with space mesh ∆a = 0.001.

Technical Details
Throughout, when BV functions are considered, we refer to a right continuous representative. We now recall the following elementary estimates on BV functions. Passing to the estimates on the total variation, introduce which is finite by (g) and (d).
Consider now the total variation estimates. Using Apply now (4.1) to obtain: completing the proof of (2.8). Concerning the stability bounds, (2.3) implies An application of Gronwall Lemma yields the desired estimate (2.11). Finally, the monotonicity property (2.12) directly follows from (2.3).
The following elementary lemma is of use below.
Lemma 4.1. Let H, K ∈ R + and assume that the numbers Proof of Theorem 2.4. Fix a time T so that which is finite by induction. Lemma 2.2 then ensures existence and uniqueness of a solution to (4.11)-(4.12) for any k > 0. By construction, the choice (4.10) ensures that u k (t, x) = u 1 (t, x) for all x > γ(t) and k ≥ 1 . Therefore, also α i t, u k (t,x i ) = α i t, u 1 (t,x i ) for all t ∈ [0, T ], for all k ≥ 1 and all i = 1, . . . , n.
Compute now Adding up all the components, and, recursively, .
Choosing now also T < 1/ n Lip(β) , the sequence u k is a Cauchy sequence and we obtain the existence of a map u * ∈ C 0 [0, T ]; L 1 (R + ; R n ) which is the limit of the sequence u k , in the sense that lim To prove that u * solves (1.1), it is sufficient to check that the boundary condition is attained. Indeed, proving that u * is a weak solution to the balance law is a standard procedure. Clearly, the initial datum is attained, since u k (0) = u o for all k.
Using Lemma 2.2, (2.3), (4.13) and (4.15), for all large k ∈ N we have b k+1 With the above choice of T , using (4.14) and Lemma 2.2, Inserting now the estimate (4.16) in the latter term above, we can apply Lemma 4.1 to the Hence, as soon as T is so small that we obtain a bound on TV b k ; [0, T ] uniform in k. This bound, thanks to Lemma 2.2, ensures that also sup t∈[0,T ] sup k∈N TV u k (t) < +∞ and, by the lower semicontinuity of the total variation with respect to the L 1 topology, also sup t∈[0,T ] TV u * (t) < +∞. Therefore, the trace lim x→0+ u * (t, x) exists for all t ∈ [0, T ].
The time T chosen above depends only on β, min ixi , on d and on g. In particular, it is independent from the initial datum. Hence, the above procedure can be iterated, extending u * to a function defined on all R + , i.e. u * ∈ C 0 R + ; L 1 (R + ; R n ) .
A repeated application of the estimate above on the intervals [(k − 1)t, kt] allows to complete the proof.