Fisheries management as a Stackelberg Evolutionary Game: Finding an evolutionarily enlightened strategy

Fish populations subject to heavy exploitation are expected to evolve over time smaller average body sizes. We introduce Stackelberg evolutionary game theory to show how fisheries management should be adjusted to mitigate the potential negative effects of such evolutionary changes. We present the game of a fisheries manager versus a fish population, where the former adjusts the harvesting rate and the net size to maximize profit, while the latter responds by evolving the size at maturation to maximize the fitness. We analyze three strategies: i) ecologically enlightened (leading to a Nash equilibrium in game-theoretic terms); ii) evolutionarily enlightened (leading to a Stackelberg equilibrium) and iii) domestication (leading to team optimum) and the corresponding outcomes for both the fisheries manager and the fish. Domestication results in the largest size for the fish and the highest profit for the manager. With the Nash approach the manager tends to adopt a high harvesting rate and a small net size that eventually leads to smaller fish. With the Stackelberg approach the manager selects a bigger net size and scales back the harvesting rate, which lead to a bigger fish size and a higher profit. Overall, our results encourage managers to take the fish evolutionary dynamics into account. Moreover, we advocate for the use of Stackelberg evolutionary game theory as a tool for providing insights into the eco-evolutionary consequences of exploiting evolving resources.

. The profit for the domestication approach. For each subplot we assume that the manager can control two variables and maximizes the profit function with respect to them for different values of the third variable. In the first subplot, the managers can control u and z and maximize the profit with respect to H; in the second, they can control H and z and maximize the profit with respect to u; in the third, they can control H and u and maximize the profit with respect to z. The profit function is concave with respect to H, u and z.

Evaluating the equilibria
The ecologically enlightened manager The ecologically enlightened manager's optimal harvesting rate is a function of the fishes' size strategy, given by ∂Q ∂H = 0.

The evolutionary enlightened manager
When u = g/(H + d), the necessary condition for the optimal harvesting rate H EV O is found by substituting g/(H + d) for u in Eq 1, taking the derivative with respect to H and setting it equal to zero.
which can be negative or positive, depending on the value of z. With our set of parameters (s = 1, g = 5, R = 4, d = 0.2, c = 5) this derivative is always negative for z ≤ 13.8. Thus, the optimal harvesting rate H EV O will always be lower than g/z − d, and, consequently, the fish size will always be equal to g/(H + d) when z ≤ 13.8. When z > 13.8, H = g/z − d, u = z and the evolutionarily enlightened approach will coincide with the one of the ecologically enlightened manager.

Stability of the equilibria
In this subsection we give a straightforward derivation of the equation governing the fish population and determine the stability of the equilibrium solutions of (S2) as a function of the parameters u and H. This is essential in order to determine whether the equilibrium solutions found for the different strategies indeed give rise to a stable fish population.
To investigate whether the different strategies will lead to feasible and stable values for the fish population, that is stable N * , we investigate the linear stability of the equilibrium N * . The adult fish population of size N t is governed by the equation where s is the fecundity rate, d the natural death rate, and R determines the effect of the surroundings on the growth of the fish and M (u) = e −du/g models the natural S 2/7 population mortality, H (u) = H for u ≥ z and H (u) = 0 for u < z. If we assume here that the individual fish chooses the same strategy as all the other fish, Eq (S2) can be derived as follows. Assume that we have two populations: a population of juveniles and an adult population of fish. The assumption that the juveniles grow linearly in time with growth rate g leads to a partial differential equation for the distribution of juveniles j(v, t) of size v at time t given as with boundary condition which accounts for suN t e −Ntu/R neonates being spawned at time t, where u is now the size of an adult fish. If we solve Eq (S3) and use T = u/g, we obtain: The adult fish now evolve according to where the contribution j(u, t) to the growth of N t is the number of juveniles promoted to adults of size u = gT and the loss is due to fishing H and natural death d. The stationary state of Eq (S2) satisfies This equation has the trivial solution N triv = 0 and a non-vanishing solution We can easily calculate the linear stability of the solutions. The linear stability of N triv is found by linearizing Eq (S2) around 0: N (t) = δN (t) + 0, which yields Substitution of δN (t) = e at , where is a small number, gives up to linear order in from which we find the condition for stability for N = 0 We remark that this condition coincides with the condition for N * to be negative. Hence as soon as a nontrivial solution for N exists, the trivial solution turns unstable. Performing a similar analysis on the N * solution leads to the following equation in terms of a This equation has only real solutions for a if a < 0, implying stability of N * to all in phase perturbations. For out of phase perturbations we consider complex solutions of S 3/7 the form a = a R + ia I . Substitution of a = a R + ia I in δN (t) = e at , with a small parameter yields the following two equations from which a I can be solved in terms of a R , (S11) By next taking the ratio of the two equations in (S10), we can find that a R is given in terms of a I as . (S12) By one more substitution, namely of Eq (S11) into Eq (S12), we arrive at . (S13) We refer to the right hand side of (S13) as f (a r ). If Eq (S13) does not have any positive solution for a R , the solution N * is stable, otherwise there will be oscillatory instabilities. Let us first remark that in order for instabilities to occur it is necessary that −1 + log suP (u)M (u)

H+d
> 1 since otherwise the first equation of Eq (S13) will never have a solution. Next we notice that a solution of Eq (S13) is expected positive when the argument of tangent function is between π/2 and π, as then the second term as a whole is positive. We can bound the conditions for a stable equilibrium by the following inequality 1 + 1 + π 2 g 2 4(H + d) 2 u 2 < log suP (u)M (u) H + d < 1 + 1 + π 2 g 2 (H + d) 2 u 2 (S14) Of course, the boundary between the stable and unstable regime is determined by setting a R = 0 in Eq (S14), which leads to (S15) In Figure S2 we plot the result of the stability calculations. More insightful is Figure S3 in which we plotted the different stability regions. From Figure S4 it is clear that for small values of H (H < 0.3) and u ≥ 30, the obtained solutions will not be stable. Numerical investigations show oscillatory behavior in this region. Eq (S2) is coupled to the equation of the yield rate Y , which is equal to P + cH, which can be controlled by In the left-hand panel we plotted the solution of Eq (S8) for two different values of u to study the stability of N = 0. The orange (dashed) curve corresponds to u = 6 and the red (solid) curve to u = 10. As the orange curve crossed the line y = x (or rather f (a r ) = a r ), the trivial solution N = 0 is stable for u = 10, but unstable for u = 6. In the right-hand panel we plotted the stability curves for the N * given by Eq (S13). Again the red (solid) curve (corresponding to u = 7) does not intersect the blue line and is therefore stable, whereas the orange curve corresponding to u = 10 does intersect the blue line f (a r ) = a r and is therefore an unstable solution. All other parameters where kept the same in both plots. More precisely z = 3, g = 3, d = 0.2, s = 1, and H = 1. Hence we find that the solutions are stable in a rather limited interval of u values u ∈ [6, 10]. Of course this also depends on the values of the other parameters.
the fisherman. The rate of biomass harvested at time t is, which also follows easily from the interpretation of the distribution of fish can be obtained as where we used that the total yield is the yield of the adult fish and the juveniles for the first equality and the expression for j(v, t) from Eq (S4) to arrive at the second equality. When we next evaluate expression (S16) in the limit that t → ∞, that is for N t → N * , with N * given by Eq (S6), we find the equilibrium value for Y as Indeed Eq (S17) leads to exactly the same profit function Q = Y − cH as the one given in Eq. 4. To summarize, we found in this subsection how taking the viewpoint of fish distribution leads to a natural derivation of the evolution equation for the fish population as well as for the expression for the profit Q. Moreover, we have found that the stability region for N * stretches over a rather large part of the H values. Only for very small values of H, instability in the fish population may arise, which could lead to modified conclusions regarding the optimal strategy. In practice this would not occur, as typical the optimal values for H are sufficiently large. The green area denotes the boundaries as given by Eq (S14). The blue curve is numerically calculated stability boundary curve following from Eq (S15). To the right of the blue curve the stationary solution is stable until the dashed magenta curve is reached, which is defined by N * = 0. When H is too large the fish will go extinct and N = 0 is the only stable solution. The brown dashed curve is the optimal strategy curve for the fish u = g/(H + d).  S4. The stability regions in the (H, u) plane for the parameter values z = 3, g = 5, s = 1, d = 0.2. The green area denotes the boundaries as given by Eq (S14). The blue curve is numerically calculated stability boundary curve following from Eq (S15). To the right of the blue curve the stationary solution is stable until the dashed magenta curve is reached, which is defined by N * = 0. When H is too large the fish will go extinct and N = 0 is the only stable solution. So the stability region for u determined by u = g/(H + d), is for H ∈ [0.1, 0.8].