A bio-economic application to the Cape Rock Lobster resource using a delay difference modelling approach

In many species, like the Cape Rock Lobster (Jasus lalandii), the life cycles of males and females differ. This may motivate the use of two-sex models in a stock-assessment analysis. It is also true for this resource, that juveniles do not reach sexual maturity immediately. Therefore a delay-difference model is appropriate. In this study we follow a bio-economic approach and use a two-sex delay-difference model to determine a maximum economic yield strategy. Thus we determine an economic optimum steady state solution at which to harvest this resource subject to the biological constraints of the species.


Introduction
In our attempts to model the population structure of species it is often necessary to include the entire life history of the population, as some events do have an immediate and homogeneous long term impact.A natural way to model these long term effects on the population dynamics is by using delay difference equations.The use of such discrete time population models is popular in fishery management strategies as is the case, for example, with the Pacific salmon fishery [5].
The standard model reads, where N t is the density of the population at time t and F (N t ) is a recurrence function.
Note that the population density at time t + 1 depends entirely on the population density at time t.
The population structure of any developed species, like birds or fish, includes both sexually mature and immature individuals, depending on their age.Thus, sexual maturity should E Roos be incorporated into the difference equation.A delay of T years in the model implies that it takes a newborn T + 1 years to reach sexual maturity.In some cases recruitment (that is the number of newborns that reach spawning age) is assumed to be independent of the current population and it is also assumed that a proportion of the current population survives to the next year.Assume also that the recruitment depends entirely on the density of the population T + 1 years earlier.Then a delay difference model for the density is where σ (0 < σ < 1) is the survival coefficient and F the recruitment function [3].
In some species one distinguishes between genders at certain stages in their lifetime.Such species are called dioecious species.Important considerations in the life cycles of sexes are mortality, growth and development, and fecundity.As discussed by Caswell and Weeks (1986) the male mortality of humans almost always exceeds female mortality.They also addressed some other significant differences among sexes in nature.Often females reach maturity later than males.For example, the mature female turtle Pseudemys scripta is 2.7 -3.2 times the mature age of the males, while in some freshwater fish species females are about 1.3 times the age of males at maturity [1].Other experiments showed that sex-specific predation rates of freshwater copepods (Diaptomus) occur [10].
Sexual dimorphism in fish species is of significance in population dynamics, as for instance the growth rate of the two sexes may differ.It has therefore become important to extend the standard delay model to a two-sex model, as was done by Cruywagen (1996).He investigated the effect of the delay time and also that of a constant effort harvesting strategy, on a population's stability.We discuss some of these results here.In §1. 1 we consider the recruitment model that was proposed by Cruywagen (1996).As is the case for the monoecious species (not shown here) periodic solutions result for specific parameter values.These bifurcations to two-, four-or eight-cyclic solutions are illustrated in §1.2.
We include harvesting in the model in §1.3 and in §1.4 a general bio-economical steady state solution is determined.
In §2 an application to the Cape Rock Lobster (Jasus lalandii ) follows.We formulate a biological model, based on the work of Cruywagen (1995), in §2.1.In §2.2 the use of the available data in the model is explained and a maximum likelihood parameter estimation procedure is used to estimate some outstanding biological parameters in §2.3.Once all biological parameters are known, it is possible to investigate a bio-economic modelling solution.Some bio-economic equilibrium results are given and discussed in §2.4.Conclusions follow in §3.

The Recruitment Model
Stock-recruitment analysis concerns the empirical relationship between the spawning stock (or mature population) and the recruitment of individuals to the spawning stock.Usually it is the spawning stock that is controlled by fishery management procedures.A thorough knowledge of all life stages prior to joining the mature population is thus necessary to be able to predict the effect that a possible change in the spawning stock might have on recruitment.A stock-recruitment curve attempts to describe how the mean recruitment rate varies with stock size.
Suppose S is the stock size, then the Ricker stock-recruitment function has the form where c denotes the initial slope of the curve and d is some constant [12].Now consider a two-sex model.Let N m t represent the number of males at time t, while N f t represents the number of females at time t.Assume the survival factors for male and females to be σ m and σ f respectively and the time delay factors u and v respectively of the male and female population.Then a non-harvested population model with a sex-structure is given by where H is a function representing the number of recruits that will develop into males and G is a function representing the number of recruits that will develop into females [7].
For dioecious populations the relative proportions between males and females become important and will be modelled here.First we wish to define the recruitment function.
Let the function ψ represent the number of zygotes (first form of newborn) produced as a proportion of the number produced when an unlimited number of males are available.
Then the number of zygotes, n(t), produced by the mature population in year t is See Bergh (1991) for a discussion of the properties of the function ψ.He proposed using the modified harmonic mean marriage function where κ is positive and real.This function is an adaptation of the reliable harmonic mean marriage function that is used to model mammal populations [4].
Here we use (4) to define the recruitment functions H and G.It is likely that in certain populations the number of male to female zygotes are not equal and this is addressed in the model.Let θ and (1 − θ) be the proportions of the zygotes that are males and females respectively.Since there is a delay of u years for males and v years for females before the zygotes become mature, the population density at each year is important and should be incorporated.In terms of the model above this implies that where F (n) is the stock recruitment function [7].
and  The dynamical solution at equilibrium of the two-sex model (3) using equations (5) and (6).(i) Male stock against the parameter c. (ii) Female stock against the parameter c, where θ = 0.6, σ = 0.0, and male and female delay factors are u = v = 0.0.

Periodic Solutions and Bifurcations
We now address the assumption made by Botsford (1992), namely that all individuals contribute equally to recruitment.We repeat the method used by Botsford, but this section goes further than his studies and differentiates between the sexes.We use ( 5) and ( 6) in the two-sex delay difference model (3) and investigate the effect of an added age structure on the model by plotting the numerical solution of N m t and N f t over a number of time steps (for example 100 steps) after it has reached equilibrium (we take it that the system has reached equilibrium after 100 steps) as function of c.We first wish to investigate the effect of the proportion of male to female zygotes on the steady state solution.Secondly, we introduce a non-zero survival rate and illustrate how the biological process changes when a fraction of the parental stock survives to the next year.Thirdly, we also investigate what effect the delay factors on recruitment will have on the solution.

Example 1
Figure 1 shows the dynamical solution of (3) for the case where u = v = 0 and k f = 0.4.Note that a series of bifurcations occur.These bifurcations are simultaneous for the male and female populations.The regions for periodic solutions are listed in Table 1, with the interval for a stable steady state being 0.0 < c < 3.47.In this example we see that a smaller proportion of zygotes develop into females (θ = 0.6), which causes a female steady state solution which is lower than that of the males.Note also that a zero steady state for the stock size occurs in the lower interval of c, which is certainly not desirable.Table 1: Approximate bifurcation positions for both male and female populations of (3), using equations ( 5) and ( 6).
(i)  The dynamical solution at equilibrium of the two-sex model (3) using equations (5) and (6).(i) Male stock against the parameter c. (ii) Female stock against the parameter c, where θ = 0.5, σ = 0.0 and the delay factors are u = v = 0.

Example 2
The proportion of zygotes that develop into females also has an effect on the stability region, as will be noted by the difference in the stability interval, as shown by Figure 2.
Here the only difference in parameter values from those in Example 1 is that θ = 0.5.The dynamical behaviour is also listed in Table 1.The interval for a stable steady state is now 0.0 < c < 3.44.Note also that the male and female populations are now equal in size, as expected.Again, a zero stock size occurs for similar values of c.

Example 3
Next we wish to illustrate the effect of the survival factor on the stability intervals.In this example we have introduced an equal survival factor of 0.05 to the male and female stock in (3).The remaining parameter values were taken as in Example 2. The interval on c for a stable equilibrium population level here is 0.0 < c < 3.46, while we find chaos in the region c > 4.13.Thus, the introduction of survival increases stability.Compare also the regions for periodic solutions of Examples 2 and 3 in Table 1.In Figure 3 this effect, of an increased survival factor on the stability, may be illustrated even better when we increase the survival factor to 0.5.The results for the female stock are identical.Here it is interesting to note that the interval where a zero stock size occurs, is in a lower interval of c, compared to the interval of Examples 1 and 2.

Examples 4 and 5
Here we use the same model and add a delay factor of u = v = 1. Figure 4(i) illustrates the stability behaviour of the system with a survival factor of 0.5.The critical value after which chaos occurs is c = 3.45.To investigate the effect the delay factor has on the biological behaviour of the stock we now increase the delay to u = v = 3.Here the critical value is c = 2.69, which indicates a smaller range of c values where stability in the behaviour exists (see Figure 4(ii)).Thus, increasing the delay factor destabilizes this equilibrium.Note also that outside the stability region chaos appears immediately, without an interval where periodic solutions occur.

Harvesting
We introduce here the model of Cruywagen (1996).Take N m t and N f t to be the preharvest densities of the breeding male and female populations at time t.Assume that the harvesting takes place during the time interval [t, t + ∆t], where ∆t is small.Assume (i)  further that reproduction takes place during the interval [t + ∆t, t + 1 − ∆t] and that recruitment occurs during the time interval [t + 1 − ∆t, t + 1].Let h m t and h f t be the male and female harvest functions, respectively, then the post harvest densities at time t are The general two-sex delay difference model (3) then becomes where H(S m , S f ) and G(S m , S f ) denote the male and female recruitment functions as before.In this model it takes males u + 1 years to reach maturity, while it takes females v + 1 years.We expect no recruitment when either the male or female population has been depleted, and also that the recruitment functions should remain non-negative.Therefore Cruywagen (1996) assumes the following properties for the functions H(S m , S f ) and G(S m , S f ): Cruywagen (1996) has shown the existence of a steady state solution and has conducted various analyses on how the delay periods, the survival factor and the constant harvesting effort, affect the population stability at the steady state.

Optimal Harvesting Policy
Finally we introduce bio-economic aspects to the model.Clark (1990) took the total value of all future net revenue, namely the present value and determined the equilibrium population for a monocious species that maximizes the present value.This maximum is subject to the biological model.
We now extend Clark's work for discrete single-sex delay models to discrete two-sex delay models.We apply this model and determine an optimal harvesting strategy, by following Clark's solution and investigating a maximum economic yield strategy by applying the standard Lagrangian optimization method.
Given a feasible harvest policy we model future stock levels using equation ( 7).Thus with initial population levels, N m 0 and N f 0 , and past survival levels S m −w and S f −w for w = 1, 2, . . ., β known, the population levels can be determined.
The sole owner's bio-economic objective is to maximize the present value function where with i the discount rate and Π k the profit in period k.Let C m and C f be the total cost functions of catching male and female fish respectively.We assume that where r is now defined as the cost per unit effort and q m and q f the catchability coefficients for the male and female populations respectively.If we take p m and p f to be the price per unit male and female fish, then the net revenue function is The optimization method with Lagrange multipliers may be used to determine the steady state population for optimal harvesting.The Lagrangian for the objective functional defined above, subject to the two-sex delay model, is where λ (1,k) and λ (2,k) are the Lagrangian multipliers.
We wish to find an equilibrium solution that satisfies these equations.Such an equilibrium is an optimal harvesting level in the management strategy.Also, at the steady states we may drop the subscript, k, on the populations sizes and hence also on Π.Thus at the steady states the population sizes satisfy the equations where as was discussed in §1.3.

E Roos
Then Use equations ( 18) and ( 19) together with the equilibrium equations ( 9) to solve for N m , N f , h m and h f .When equation ( 8) is used we find that Similar expressions are determined for the females.Equations ( 18) and ( 19) are now: 2 Application to the Cape Rock Lobster In this section the two-sex delay model is applied to the Cape Rock Lobster Resource (Jasus lalandii ).

Biological Model
Assume the total number of zygotes in year τ to be Z(S m τ , S f τ ), where S m τ and S f τ is the post harvest male and female densities.Assume also that a proportion, θ, of these zygotes will develop into males.It is observed in nature that female juveniles take longer than males to reach sexual maturity.We therefore assume that v > u.The general two-sex model (7), where In §1.1 the adapted harmonic mean marriage function was discussed along with the Ricker stock-recruitment functions.Here we propose an adapted version of the Ricker-type function to describe recruitment.For the rock lobster, we assume that enough males are available so that where κ, β, δ and ω are positive constants.In general, with this proposed function, recruitment is higher than that in the original Ricker function in §1.1 and decreases asymptotically to κ instead of zero.In order to have a zero recruitment with zero stock it is necessary to shift the entire function to the right with the introduction of δ.Also, as with the Ricker function, with an increase in stock size the mortality rate increases for high egg productions.
As the available data are in terms of biomass, rather than in terms of population size, we have to adapt the above model ( 22) accordingly.We use the method followed by Cruywagen (1995) and perform this conversion using the Ford-Walford model [8].That is, assume that the growth in mean body weight at age, a, may be described by the linear relationship where w m (a) is the weight at age a, and α m and ρ m are empirical constants [8].The equation for females are analogous, with f as the subscript instead of m.
Assume M (t) to be the total biomass of males at maturity (age u + 1 and older) and similarly F (t) to be the total female biomass at age v + 1 and older.Following Cruywagen's (1995) approach, it is now possible to express biomass, in terms of numbers multiplied by body weight, namely S m t+1 (a)w m (a), where S m t+1 (a) refers to the number of males aged a, in the mature male population in year t + 1, while S f t+1 (a) is the equivalent for females.If h m t+1 (a) represents the number of males of age a that have been harvested, and m t+1 (u + 1) represents the number of the previous year's recruits that have been harvested, then E Roos Similar relationships may be obtained for females.After substituting (24) in ( 23) we have , where and similarly for females.
Average weights, w m (t) and w f (t), for mature males and females respectively in year t, are defined as This is similar to what was done by Cruywagen (1995).The equations then reduce to At the population's pristine level, we assume and the population levels in numbers are If we ignore the catches (because we are at the pristine state) can use the numbers model ( 21) and biomass model (25, 26) to derive the relationship between the average male weight at the pristine state, w m , and the natural male survivorship parameter, σ m , namely For the females the equation is Equations ( 25) and ( 26) at equilibrium may be used to determine the relationship between the male and female carrying capacities, [6].By using the last three equations it is possible to reduce the number of parameters to be estimated, as three of them are given here in terms of the others.

Data
As this is a two sex model we seek data sets on both males and females.We now assume that the biological parameters u, v, σ f , σ m and θ are known.
The values of the delay times are taken as u = 6 and v = 7 (according to the base case in the study by Cruywagen (1995)) and there is no reason not to take θ = 0.5.
Johnston and Berg (1993) were able to estimate the male survivorship parameter to be in the range 0.82 to 0.95.Here, the values of σ m = 0.90 and σ f = 0.92 are adopted, as were used by Cruywagen (1995).
The parameter K F is expressed in terms of the others, by using equation ( 29).
Other data sets required are the catch per year, catch per unit effort, Ford-Walford parameters and average weight per year.Parameter values in the model are the male capacity, K M and the recruitment modelling parameters κ, δ, ω and β.
Due to an agreement with Marine and Coastal management of South Africa we do not publish specific data here, but rather discuss the process that was followed.Data that are available for certain years are the total catch per season, the percentage males in the catch and the catch per unit effort series.Also available are tagging data of different fishing areas during different years.
The total catch data is available for the period 1900-1995, together with the proportion of males in the catch.(We assume the population to have been at its pristine level prior to 1900.)Refer to Cruywagen (1995) for a detailed discussion on the catch sex ratio.Note that it is common to find a very high percentage of the catch to be males.It is also necessary to note that catches dropped drastically after the very high catch rates during the 1950's and 1960's.
The catch-per-unit-effort (CPUE) data are available since 1975.Note here the catch per unit effort data for males and females are in the same proportion as in the catch.Until the 1991 season a legal minimum carapace length of 89 mm was enforced.In the 1992 season the legal minimum carapace length was dropped from 89 mm to 80 mm.Since the 1993 season it has been dropped even further to 75 mm, mainly due to the unavailability of larger lobsters and the market preference for smaller lobsters.As a result, with the same effort as before, larger catches were taken which thus lead to an increased CPUE.
We use the Ford-Walford growth parameters as was calculated using the tagging data [6].
Next it is also necessary to determine the average weights of males and females per year, starting at the pristine state.Cruywagen (1995) also reported that a preferred method of obtaining size-frequencies is by using the results of diving experiments, rather than the examination of catches, as these are selective catching.Such experimental results are reported in [6].
We follow Cruywagen (1995) and assume that the decline in the average weights is proportional to the total of catches over years.He uses a linear interpolation method that reduces the weights w m and w f (1900) towards the average weights of males, w m (t) and females, w f (t) per year.This decline is proportional to the total catch taken until E Roos that year, namely, , where only the pristine average weights, w m (1901) and w f (1901) and the 1995 average weights, w m (1995) and w f (1995), are known.Equations ( 27) and (28) were used to calculate the pristine average male weights.In terms of lengths that is 101 mm for males, and 71 mm for females.Take the 1995 average mature male weight to be that of a 75 mm animal and that of the females of a 68 mm animal.(These lengths agree well with the experimental data for males, but differ to some extent from those of the females.)The calculated average weights for males and females are given in Figure 5.It is noticeable that the male curve has a very steep decline compared to that of the female curve.The decline is the direct result of very high catches of earlier years and the fact that most of the catches landed, consists of males as they are naturally larger (in size) than the females.The parameter estimation procedure that is described in the following section was used to estimate the other six parameters, K M , K f , κ, δ, ω and β.Equation ( 29) is used to determine K f .Three parameters, κ, δ and ω were fixed at some chosen values and by fitting the model the two remaining parameters could be estimated, namely β and the male carrying capacity K M .The Powell optimization method was used in the parameter estimation process [11].Table 2: The results of the model application on the Cape Rock Lobster with various parameter choices in the model.The parameters κ, δ and ω are chosen while the algorithm searches for K m and β.

Estimation of Biological Parameters
We use the maximum likelihood parameter estimation procedure to estimate the male CPUE.Assume that the CPUE is proportional to the average between the pre-and post harvest male population sizes, namely where q m is the male catchability coefficient and E(t) the effort level at time t.Similarly the female CPUE is modelled by where q f is the female catchability coefficient.
It is standard practice to assume that the observed catch-per-unit-efforts, CPUE m (t) and CPUE f (t) have a lognormal distribution [8].The model predicted CPUE may then be as CPUE m (t) = CPUE m (t)e ε , with ε = N(0, σ 2 ) a normal distribution with mean 0 and standard deviation σ.A similar expression holds for the female CPUE.
By fitting modelled CPUE series to the observed series the unknown parameters in the delay model (25, 26) could be estimated.

Define the likelihood function as
where p refers to the number of years for which data are available.We have discussed earlier why the catchability coefficient differs before 1992 and thereafter.Thus take q m,1 to be the catchability coefficient for males during the period 1990-1991 and q m,2 to be the coefficient for the period since 1993.Similar notation is used for the females.

E Roos
At the maximum for the log-likelihood function the catchability coefficients satisfy with standard deviations The Powell minimization method was used to perform the CPUE fit to the observed data.Various choices for the three parameters κ, δ and ω of the recruitment function were investigated and compared (see Table 2).As it is not easy to estimate all unknown parameters by the fitting procedure we used a trial and error process to decide on realistic parameters for the recruitment function in (22).A typical set is κ = 3.0×10 7 , δ = 1.0×10 6  and ω = 5.0 × 10 −8 .The best fit for this set estimates β = 610 and K m = 4.614 × 10 8 .The recruitment curve with these parameter values is given in Figure 6. Figure 7 shows the predicted male and female stock size of the model ( 25) and ( 26) for each year.Note the general decline since 1900 in stock size for both the males and females.Again, as was the case with the average weights predicted, the difference in the male and female curves is due to the difference in male and female catches.

Bio-economic Modelling
It is possible to determine expressions for the optimal bio-economic population and catch levels in terms of the production function given in §2.1, namely equations ( 21) and ( 22).
At the bio-economic equilibrium the biological model ( 21) and equation (20) have to be satisfied.Here, and the set of equations becomes Thus we have four equations in the four unknowns N m , N f , S m , S f , as all other parameters have been estimated in the previous section.
Let us now consider a steady state solution for equations (30), using the estimated biological parameters of the previous paragraph.Take Γ = exp(−ωS f ) in equation ( 22).Then Γ → 0 when S f → ∞.To assume Γ ≈ 0 is not general, but possible for large female stock sizes.Under conditions where the function ∂Z ∂S f → 0 the set of equations (30) simplifies to the point where they may be solved analytically.The equations reduce to a pair of second order polynomials and the steady state solution thus tends to the roots of these polynomials as S f → ∞.Without describing the analytical expressions here, we rather show some results.Take the biological parameter values as in the previous section, then Table 3 shows some of these solutions for different economic values.The results of the stock sizes are given in numbers.This table gives the equilibrium stock level to be maintained in order to maximize the total present value for different economical factors.   2 (last row).N m is the male pre-harvest stock size in numbers, N f is the female pre-harvest stock size in numbers, S m is the male post-harvest stock size in numbers, S f is the female post-harvest stock size in numbers, h m is the annual male harvest, h f is the annual female harvest and Π is the annual net revenue from such a strategy.
Table 3 shows that the female stock size decreases as the price per female increases, due to a higher catch demand.On the other hand, the female stock increases with an increase in the cost per unit effort as the catch will probably decrease.As expected, the female stock size decreases as the interest rate increases, because future stock will be worth less under higher rates.
Another interesting result from Table 3 is that the male catch size is independent of the price per female fish and vice versa.

Conclusion
For species like the South African West Coast rock lobster there is a delay before newborns reach sexual maturity.This phenomenon has been addressed by using a delay difference model.For many species the population dynamics of the male and female population are different and it may become necessary to use a two-sex delay difference model in order to distinguish between the sexes, as was done in this paper.
The overlap of generations has been incorporated in the model by the use of a survival factor.As an age structure (or increase in survival) is added to the Ricker model the critical value, where bifurcation occurs, increases such that increased survival increases stability.Contrary to this we saw that with the addition of the delay factor, the critical value decreases and delay is thus destabilizing.
We used the theoretical framework of Cruywagen (1996) and added bio-economic features to the model.By applying optimal control theory we arrived at the four equations in ( 9) and (20) to solve for the four unknowns, N m , N f , S m , S f , at equilibrium.
In §2 the proposed two-sex delay difference model was applied to the Cape Rock Lobster resource.Data on this fishery are available only in biomass.Consequently a conversion from length to weight was necessary in the model.The pristine average male and female lengths were calculated from equations ( 27) and (28), and the female lengths differ slightly from the range on average length from experimental data.This difference may indicate that the model does not properly describe the population, but that the experimental data are only for certain exploited areas and may therefore only be taken as an indication.
A best fit of the modelled CPUE to the observed CPUE was performed.Shortcomings of this attempt include the drawback that the total outcome relies only on a CPUE series beyond 1975, when catches were already low.Earlier data are not available.Better model fits are anticipated in situations where a longer time series is available.
Nevertheless, the model provides a much more detailed description of the resource than the usual surplus production model.Although we have made assumptions on the recruitment function and harvesting outcome, the above model and results could be usefull in planning specific future management strategies for male and female populations.Also, environmental effects on the current stock will influence the recruitment T + 1 years in the future as the recruitment depends on the delay factors.

Figure 2 :
Figure 2:The dynamical solution at equilibrium of the two-sex model (3) using equations(5)

Figure 5 :
Figure 5: Calculated average male and female weight time series for the Cape Rock Lobster resource given in kg.

Figure 7 :
Figure 7: Predicted male and female population size in biomass for the Cape Rock Lobster.

Figures 8 (
Figures 8 (a) and (b) illustrate the fits of the modelled CPUE to the observed male and female CPUE values.

Figure 8 :
Figure 8: (a) Male modelled CPUE fit to observed data for the Cape Rock Lobster.(b) Female modelled CPUE fit to observed data for the Cape Rock Lobster.

Table 3 :
The bio-economic equilibrium results for different economic values when using the biological values given in Table