Richards growth model and viability indicators for populations subject to interventions

In this work we study the problem of modeling identification of a population employing a discrete dynamic model based on the Richards growth model. The population is subjected to interventions due to consumption, such as hunting or farming animals. The model identification allows us to estimate the probability or the average time for a population number to reach a certain level. The parameter inference for these models are obtained with the use of the likelihood profile technique as developed in this paper. The identification method here developed can be applied to evaluate the productivity of animal husbandry or to evaluate the risk of extinction of autochthon populations. It is applied to data of the Brazilian beef cattle herd population, and the the population number to reach a certain goal level is investigated.


INTRODUCTION
There are multiple reasons for modeling the growth of a population of a certain species.A growth model of a population is an important tool for understanding how environmental uncertainties affect growth, and it provides the foundations to growth control policies, slaughter control, etc.These models help to forecast populational behavior and to evaluate the risk of extinction or, in the other extreme, population explosions.From the stochastic viewpoint, they allow calculations of the probability of extinction, explosion and to achieve a particular population goal or the expected time for the occurrence of these events.These 1108 SELENE LOIBEL, MARINHO G. ANDRADE, JOÃO B.R. DO VAL and ALFREDO R. DE FREITAS populational characteristics are known as population viability indicators.The population viability analysis (PVA) is one of the central tools for conservation planning and evaluation of management options.Compared to other methods, PVA has several advantages (Boyce 1992, Omlin and Reichert 1999, Akçakaya and Sjögren-Gulve 2000).PVA provides a rigorous methodology that can use different types of data a way to incorporate uncertainties, natural variabilities and products or predictions that are relevant to conservation and management goals.The disadvantages of PVA include its single-species focus and requirements for data that may not be available for many species.
There are many works dealing with the population growth of a single population with growth rates depending on its size (see Ludwig 1996, 1999, Allen and Allen 2003, Matis and Kiffe 2004).In (Ludwig 1996(Ludwig , 1999)), the Gompertz and the logistic model for single species are considered to calculate extinction probability.In (Matis and Kiffe 2004), a stochastic (Markovian) logistic population growth model with immigration and multiple births is derived, and the equilibrium distribution of the population is approximated to provide the extinction time of a species.The use of some of these models to evaluate the economic impact of population growth control can be found in (Yamauchi 2000).In (Matis et al. 2007), it was developed a generalization of the deterministic and stochastic population size model based on powerlaw kinetics for the black-margined pecan aphid.Deriso et al. (2008) present a framework for evaluating the cause of fishery declines by integrating covariates into a fisheries stock assessment model.In (Polansky et al. 2008), it was used simulated data to quantify the relative effects of sample size, population disturbances and environmental stochasticity on statistical inference.They focused on the key parameters that inform population dynamics predictions in a generalized Beverton-Holt model.In (Jiao et al. 2009), a hierarchical Bayesian approach for population dynamics modelling was developed to simulate variability in the population growth rates of three fish species.
In this paper, PVA is used to estimate the limiting probability of a population to reach a certain level, or the expected time for the occurrence of this event.It considers the population dynamics described by the Richards model, (Richards 1959), also known as generalized logistic model.As implied, it generalizes the logistic and also the Gompertz growth model (see Loibel et al. 2006).The Richards model was first used for the growth of animals (weight × age) and only later it was applied to model the populational growth (Fitzhugh 1974).A variation of this model, proposed by Gilpin and Ayala (1973), is known as theta-logistic, and it is used by Saether et al. (2000Saether et al. ( , 2002a, b) , b) for the study of different populations.
Most of the works found in the literature relies on the logistic or the Gompertz models due to the difficulties to estimate the form parameter of the Richards model.The models with two parameters have a growth rate per capita that is symmetric near the inflexion point at mid level of the carrying capacity.However, this is not the case of populations of large mammals (see Fowler 1981).A three parameter model, such as Richards', is able to fit better to an asymmetric growth rate per capita due to the extra form parameter (see Sibly et al. 2005).In addition to the usual Richards model, we introduce an extension to take into account exogenous interventions on the population.Modeling with interventions becomes very important when we need to identify the growth parameters of a population subjected to human intervention with economic and/or ecological motivations.
In this work, we present an analysis of the Brazilian beef cattle populational growth, here considered as a single gross population with growth rate depending on the population.Brazil currently occupies the An Acad Bras Cienc (2010) 82 (4) second place in bovine herd size and in production of beef meat, only surpassed by the USA.The Brazilian bovine population has presented growth in recent years, having about 169 million head of cattle, and an increase of more than 10 million just in the last ten years (see ANUALPEC 2008).
The Brazilian cattle production in the last decades presented constant growth rates in production terms, export and consumption.Brazil has a potential domestic market for food consumption mainly for the beef meat.The domestic cattle production has always been characterized by an extensive system, since the cattle are mostly bred as unconfined herds, being relying on pastures.In recent years, with the incorporation of new technologies that aim at increasing productivity, the intensive systems of production increased in some regions, being called confined or semi-confined.However, the physical space and other production factors needed for further growth are bound to be scarce in the future, creating a saturation of the growth rate that are well described by the class of growth models here studied.
The development of cattle raising becomes even more attractive when the social-economic potential of the meat production chain is considered.The agro-industrial system of animal meat has the largest impact on the Gross National Product, and the segment denominated "steer" alone is responsible for more than US$ 30 million a year.According to Pineda (2002), a goal for Brazil to take advantage of its full potential and to definitely hold the position as one of the international main producers and exporters of bovine meat would be to keep a bovine population stable with more than 200 million animals.
Several factors contribute to the growth of bovine populations, which assign an stochastic behavior to the growth process.As a result, the use of methodologies that provide accurate estimates of bovine population growth in the country is important, as they contribute to the organization of this sector giving support to policies on internal consumption, exports and marketing, among others.In this work, we apply specifically the Richards model with exogenous interventions on the population to make estimates of the probability of the population beef cattle to achieve a particular goal.The model also allows the average time up to the occurrence of these events by taking into account the interventions on the population.
The identification problem of the Richards model without intervention and a detailed analysis of the inference of the model parameters were previously addressed in Loibel et al. (2006) by using bootstrap techniques.As in this case, the inference of the present model is treated with the use of the likelihood profile function technique that, roughly speaking, fix some parameters (the nuisance parameters) to estimate others see also (Seber and Wild 1989).Here the intervention is taken into account as part of the population that is removed, which may be brought about by fishing, hunting or slaughter in the case of animal husbandry.In these cases, we are interested in estimating the parameters of the model, not only the ones related to the free trajectory of the process {N t , t ≥ 0}, but also those associated with the intervention {h t (N t ), t ≥ 0} modeled as a stochastic process.
This work is organized in the following way: in section 2, we present Richards' populational growth model with intervention.In section 3, we present the procedure for identification of the parameters of this model employing the maximum likelihood profile function technique.In section 4, we present a Monte Carlo simulation approach to evaluate the indicators of the populational viability.In section 5, we present the results of this analysis for the data gathered over a period of 25 years (1983 to 2008) regarding the Brazilian bovine population.We calculate the probability of the Brazilian herd reaching the goal of 200 million head proposed by Pineda (2002), and the time expected for the Brazilian herd to reach this goal if an adequate slaughter policy is established.Apart from its own importance, this case study offers signposts for similar populational studies.

RICHARDS MODEL WITH INTERVENTION
We consider the size of the population as a stochastic process {N t , t ∈ N}, where N is the set of non-negative integers {0, 1, 2, . ..}.This population consists of adults and young, but only adults can be recruited to slaughter.Denoting the number of slaughtered individuals by {H t , t ∈ N}, the number of individuals who remain in the population is R t = N t − H t , where R t represents the stock of animals available for slaughter in the next instant t + 1.
To write a growth model with an intervention function (or slaughter), we consider the population size N t and N t+1 in two consecutive instants t and t + 1, respectively, and the number of slaughtered individuals H t+1 .The total number of individuals N t+1 can be written as a function of the population size in the instant t (we denote this function by f (N t )) subtracted from individuals who will be slaughtered in the instant t + 1.Then, the dynamic model can be written as We model the interventions on this population as a function {h t (R t ), t ∈ N}, where 0 ≤ h t (R t ) ≤ 1 represents the fraction of R t of a "removed" population in the instant t + 1.Then, the number of slaughtered individuals H t+1 is given by H t+1 = h t (R t )R t .Then, the dynamic model can be written as a cause-effect relation as (1) Supposing that the total population N t+1 and N t are related by the Richards growth function, then we can write that where t is an stochastic process in which each value is independent and identically distributed (i.i.d) with distribution N (0, τ −1 ) for τ > 0. It represents the random effects of environment on the growth rate of the Richards model.Parameters ρ is the intrinsic growth rate, K is the carrying capacity, and q = 0 is a form parameter that changes the shape of the growth curve.
Replacing the equation ( 2) in (1), we can write the complete model in the presence of interventions as: To elaborate the intervention function, note that it should depend on the available population, and an admissible model for h t (R t ) is a linear piecewise model given by: where R 0 > 0 is the minimum viable population size, below which the population is considered extinct (call it "realistic" extinct).R 0 = 2 is certainly a possible choice because, with less than two individuals, there is no reproduction, but we can choose larger values for R 0 , if there are Allee effects.The rate 0 < h 0 ≤ 1, is the maximum intervention rate due to the technological and physical limitations for the intervention action.The equation ( 4) can be written as: where t is a non-correlated i.i.d.stochastic process with distribution N (0, η −1 ), with η > 0. We assume that t and t are independent processes; therefore, we have that E( t ω t ) = 0.In (5) we have Considering the model in ( 5) for the rate of intervention h t (R t ) as a function of the available population R t , we have that the conditional distribution density The equation ( 3) can be written as: We denote , and also set y h,t = q ln(N t+1 /N t ) and z t = q t .We can, then, rewrite the equation ( 6) as a linear relation: with α 0 = ρ and α 1 = −ρ/K q .From the fact that t is a i.i.d stochastic process with distribution N (0, τ −1 ), it follows that z t is also i.i.d stochastic process with Gaussian density with zero mean and variance τ −1 0 = q 2 τ −1 .We, then, conclude from the equation ( 7) that the conditional density P(y and variance τ −1 0 = q 2 τ −1 .We use this characterization to formulate the identification problem for the Richards model.

MODEL FITTING
In order to formulate the problem of fitting the Richards model with intervention, by maximum likelihood method, it is necessary to know the joint distribution of the processes {N t , R t , h t (R t ), t ≥ 0}.If we know the stochastic processes R t and h t (R t ), we can recover the stochastic process N t using the relation We begin the evaluations with the joint distribution of the processes {y we have, therefore: To write the joint density for the original processes, P(N t+1 , h t (R t )|N t , R t ), we should consider the Jacobian of the transformation y (q) h,t = q ln (N t+1 /N t ) given by: The joint conditional density for {N t+1 , h t (R t )} may be calculated, since h t is a function only of R t , by: Considering the model in ( 5) and ( 7), we can write the joint conditional density (10) for {N t+1 , h t (R t )} as: With ( 11) we can construct the likelihood function used to estimate the parameters of the model.

LIKELIHOOD FUNCTION
Adopting the matrix notation , we conclude that the likelihood function for the Richards model with stochastic intervention is given by: and In the equation ( 13), we are using the notation Y (q) = [y h,1 , . . ., y n ] and X (q) = [1, N (q) ], V = [1, R], with 1 being a n × 1 vector of ones.Due to the separability of L(α, φ, τ 0 , η, q|D) regarding the parameters (α, τ 0 ) and (φ, η), we can write the log-likelihood function as: where: To make inference about model parameters it is necessary to calculate the Fisher information matrix.We adopted the matrix notation θ = (α , τ 0 ) and γ = (φ , η).Then, for a fixed q, the Fisher information matrix for θ and γ has the form: The fact that I (q) is a block diagonal matrix implies in the independence between the estimators θ and γ .We use the likelihood function for the Richards model ( 15), for a fixed value of the parameter q to calculate the maximum likelihood estimates (MLE) of the parameters θ (q) = (α (q) , τ (q) 0 ), conditioned to a particular value of q.We denote these estimators by θ (q) = ( α (q) , τ (q) 0 ).Substituting these estimators in equation ( 15), we have the profiled likelihood function for parameter q given by The equation ( 18) is a modified version of the profile likelihood function for the Richards model, Loibel et al. (2006), which includes the presence of the intervention in the vector We should find the maximum profile likelihood estimate (MPLE) q * of the parameter q evaluating the function ( 18) for many values of the parameter q ∈ [q min , q max ], and finding q * that maximizes the function (18).
ESTIMATION OF THE PARAMETER γ (q)   The estimators of maximum likelihood for φ and η are calculated by maximizing the function l(γ |D) given in ( 16).Therefore, we have: The variance for φ can be calculated straightforward by V ar ( φ) = η −1 (V V) −1 .The relation R 0 = −φ 0 /φ 1 and h 0 = φ 1 (K + φ 0 /φ 1 ) are used to estimate the parameters present in the slaughter function h t .The variance for the vector of original parameters ψ = (h 0 , R 0 ) can be estimated with the observed information matrix I (φ) = η(V V) calculated with the observed information matrix I (φ) and the Jacobian of the transformations R 0 = −φ 0 /φ 1 and h 0 = φ 1 (K + φ 0 /φ 1 ) (to details, see Seber and Wild (1989)), given by The Fisher information matrix for the original vector of parameters is given by the variance for the estimates of the original parameters ψ = (h 0 , R 0 ) is given by and the variance for η is given by V ar ( η −1 ) = 2/(n − 1) η 2 .

ESTIMATION OF THE PARAMETER q
It is often of interest to test if the parameter q of the Richards model is equal to a known value q (0) , that is, H 0 : q = q (0) .In this way, we can easily obtain from the profile a likelihood ratio (LR) statistics w = 2 {l P (q * |D) − l P (q (0) |D)} to test q = q (0) , which has an asymptotical Chi-square χ 2 1 distribution with one degree of freedom, (Barndorff-Nielsen and Cox 1994).From this result, we can produce a large sample confidence interval for q by inverting the LR test.An approximate confidence interval for q is, then, readily obtained from {q | l P (q|D) > l P (q * |D) − 1 2 χ 2 1 (δ)}, and the accuracy of this approximation comes from the fact that Pr w ≥ χ 2 1 (δ) = δ + O(n −1/2 ).

CALCULATION OF VIABILITY INDICATORS
In general, the growth models are adjusted to allow the calculation of viability indicators of a population, such as the probability of the population to reach a determined upper level defined by G = {x ∈ R + : x ≤ g 0 }, 0 < g 0 < ∞, or an upper level set G = {x ∈ R + : x ≥ g 1 }, 0 < g 1 < ∞.These intervals may be associated to a very low level, which indicates a state of near extinction, or an upper goal, which indicates a demographic explosion or an economically and/or ecologically desired objective.

MONTE CARLO SIMULATION APPROACH
We consider the populational growth model given by equations ( 3) and ( 5), remembering that t and ω t are independent processes, with t , i.i.d.N (0, τ −1 ) and ω t , i.i.d.N (0, η −1 ).With these models, we can simulate the behavior of a population over a long period of time and generate possible trajectories for the processes {N t (R t ), j = 1, . . .J } for values of time beginning with t = T 0 , for a certain initial instant T 0 of interest.In this way, we estimate the risk or probability of the population reaching a goal by counting how many trajectories were necessary for this goal to be achieved among all of the J generated trajectories.We denote the time elapsed to achieve this goal by T G , with T G = inf{k > T 0 : N k ∈ G|N T 0 } and the number associated to these occurrences by #{T G < ∞}.Therefore, we have the probability estimate of the population reaching a certain goal by: In some cases, we may be interested in calculating not only the probability of the population to reach a certain state, but also to evaluate the expected time for the population to reach a certain state.The time necessary for the population to reach a goal is a stop time T G and its expected value can be estimated using the generated trajectories by simply evaluating the time to reach the level to each realization T ( j) G , j = 1, 2, . . .J .We, then, have the estimate: We should adopt a very large number of J trajectories in view of the law of large numbers to assure the convergence of the estimators (28) and (29) for the true values.Another difficulty found in the practical application of this technique is that, depending on the intervention policy adopted {h t (R t ), t ≥ T 0 }, the goal may become unreachable or perhaps reachable with a very low probability.It is necessary, therefore, to establish a limit for the maximum simulation time that will be associated with infinity.We denote the infinite time by T ∞ , that is, {T k < ∞} ≡ {T k < T ∞ }.Then, in the simulation process N t , we have N t < G for t ≤ T ∞ , and we consider that the goal was not reached.

DETERMINISTIC INTERVENTION
There are situations in which we can consider the intervention h t as a deterministic function.One of these situations occurs when we want to evaluate the impact of an intervention, for example, in an event where a government policy establishes that the removal of a population should be done in a certain period of the year.An estimate of the mean removal can be established, and we wish to evaluate the impact of this average removal on the population.Another situation occurs in long-term studies when we want to evaluate the impact on the population when a mean removal is carried out over a long period of time.
In these cases, the estimates of the parameters and the risk of extinction are obtained for a predetermined intervention policy.Assuming that {h t (R t ), t ≥ 0}, which means a deterministic function of R t that represents the proportion of "removed" animals at time t + 1, we can write the Richards model with interventions in the same way as in (3), and estimates of the parameters θ and q can be calculated as in the section for stochastic interventions.When the intervention h t is deterministic, the recursive equation ( 2) can generally be written as The subscript h indicates that, somehow, the function f h (•) depends on the action h.In the model proposed in this paper, the intervention h affects the maximum likelihood estimators of the parameters of the function f h (•), which are calculated using the equations ( 19)-( 21).The process { t , t = 0, 1, . ..} is independent and identically distributed, with E( t ) = 0 and E( t t+s ) = 0, if s = 0 and E( 2t ) = τ −1 .We can adopt a discrete state space model for the problem of population growth whose continuous states space E can be decomposed into the disjoint union The probability of the population to reach a certain state I j at instant t + 1, given the records of the population N t = {N 1 , N 1 , . . ., N t } and of the intervention H t = {h 0 , h 1 , . . ., h t }, is written as follows: where From (30) we see that the random process {N 0 , N 1 , . . ., N t } is function of the basic random process {N 0 , 1 , . . ., t−1 }, which is independent of t .Hence, t is also independent of {N 0 , N 1 , . . ., N t }.Under these conditions, we can write that If N t ∈ I i , for some I i ⊂ E, we can adopt a Markov chain model for the problem of population growth with the probability of transition between states I i and I j , p i j = P(N t+1 ∈ I j |N t ∈ I i ), given by: Let P be a square matrix with the elements p i j , defined for all {(i, j), 0 ≤ i, j ≤ N E }.Then, P is called Markov matrix over E if it satisfies the following properties: We are interested in calculating the probability of the chain leaving a state I i at instant t and reaching a state I j at instant t + s, s ≥ 0. Therefore, we should calculate the transition probability matrix in s steps denoted by P s .For any t ∈ [0, T ], the probability of the chain leaving the state I i for state I j in s steps is P(N t+s ∈ I j |N t ∈ I i ) = P s i j , for all I i , I j ⊂ E and s > 0. To calculate this probability, we can use the Chapman-Kolmogorov equation.Cienc (2010) 82 (4) In many problems, we may be interested in calculating the probability of a population reaching any state above I j .Thus, denoting M = N E k= j I k , we have:

An Acad Bras
In other words, the probability of reaching M is the sum of columns l, so that l ≥ j on line i of the matrix P s .
For the calculation of the transition matrix using the Richards model with intervention, we denote F h (N t ) = (ρ/q) (1 − (N t /K ) q ) and rewrite the model (3) as: where If I i and I j are any two states of space of states E, and assuming that I j = [ A j , B j ], then an interval [a j , b j ] exists so that we can calculate the transition matrix in one step with elements p i j = P(N t+1 ∈ I j |N t ∈ I i ) and rewrite the equation ( 34) by: to find [a j , b j ].We can use the equation (37) to write: Then, we can write that With equations ( 38) and (39), we can write that Assuming that t is a process with Gaussian distribution N (0, τ −1 ), we can calculate p i j from the standard Gaussian accumulated distribution function with: This form of calculating the matrix of transition probabilities is a discrete approximation, since we consider the probabilities of transition between the intervals I i and I j independent of the values that the population may assume within the intervals.The accuracy of this approach depends on the discretization refinement of the space of states.
Given that N 0 ∈ I i , we wish to calculate the expected value E(T ), so that N s makes the first passage for the state I j .In practice, there are many situations where the calculation of first-passage time is important.In the case of control of populations that grow explosively, we want to calculate the time expected for the population to reach a point that justifies intervention, such as controlled slaughter.
DEFINITION 1. First-passage time: Let I i and I j be two states of a Markov chain {N t , t ∈ I + } with space of states E. If, for a given time t, N t ∈ I i , the first-passage time of the state I i to state I j is defined as: DEFINITION 2. For any s ∈ I + , the probability function of T i j is defined by: By convention, we adopt g i j (0) = 0.
DEFINITION 3. Expected value of first-passage time: PROPOSITION 1.For every pair of states (I i , I j ) ⊂ E, we have: PROOF.Using Definition 2 for t = 1, we have: Therefore, Since P(T i j = t) = P(N t ∈ I j |N 0 ∈ I i ) and considering all of the possible states I k with k = j, we have: We, then, have: An Acad Bras Cienc (2010) 82 (4) PROPOSITION 2. E = N I k=1 I k with I i I j = ∅, ∀i = j; let N 0 ∈ I i and μ = (μ i j , 0 ≤ i ≤ N I e i = j) be the expected times to reach state I j , starting from any I i state.We can calculate μ by resolving the system of equations given by: where I is the matrix identity (N I − 1) × (N I − 1), 1 = (1, . . ., 1) is a vector (N I − 1) × 1, and Q is the matrix ((N I − 1) × (N I − 1) obtained by removing the line and the column j from the matrix P.
PROOF.We will consider that: We can write that: Replacing ( 54) in ( 55), we have: Establishing a j and writing (57) for μ = (μ i j , 0 ≤ i ≤ N I e i = j), we have the system of equations given by: The Markov chain modeling can only be used together with a constant slaughter rate, and it is more appropriate for evaluation of scenarios in long-term planning.On the other hand, the Monte Carlo simulation approach can handle more general slaughter functions but it also requires a larger computational effort.A comparison between these two approaches, for constant slaugher rates is presented in the Results section.

RESULTS
In this study, we propose an evaluation of the growth of meat production in Brazil based on data regarding the number of cattle and annual slaughter between 1983and 2008(ANUALPEC 2008)).The population and slaughter data referring to this period are presented on the graphs shown in Figures 1 and 2, in millions of heads, together with the adjusted values, according to the equations (3) and ( 5), respectively.
Table I presents the maximum likelihood estimates (MLE) and the standard deviation (SD) of the parameters of the Richards model adjusted to the population data in Figure 1, and of the parameters of the slaughter function adjusted to the data in Figure 2, as well as the 95% confidence interval for q obtained from {q | l P (q|N ) > l P (q * |N ) − 1 2 χ 2 1 }, with χ 2 1 = 3.8414.The relation N 0 = −φ 0 /φ 1 and h 0 = φ 1 (K + φ 0 /φ 1 ) was used to estimate the parameters present in the slaughter function h t .An important observation refers to the estimated values of K and h 0 that represent the carrying capacity and maximum slaughtering allowed, respectively.Note that the carrying capacity estimate K = 217.5701exceeded the maximum population level N 26 = 176.1144 million head of cattle observed in 2004 (see Figure 1) and h 0 = 0.3294 that represents an estimate of the maximum slaughter rate that is exceeding the maximum rate, which was about 0.2772, observed in 2006 (see Figure 2).K and h 0 are the two main parameters that affect the population growth.In this work, we perform simulations of population growth for different values of K and h 0 , close to estimates values, in order to estimate the probability and the expected time of the population to reach a given goal.
The graph in Figure 3 presents some of the 100 thousand simulations of trajectories of the process that were generated to calculate the probability of the population reaching the goal of 200 million heads.
In Table II, we present the probability p (Percent) and the expected time T (year) for this event to occur, considering scenarios with K increased from the estimated value to the value of 250 million heads, and h 0 increased from the actual observed value 0.21 to the estimated value 0.3294 given in Table II.
Note that the probability that the Brazilian herd will reach the goal of 200 million head proposed by Pineda ( 2002) when h 0 = 0.3294 and K = 217 million heads is only p = 0.02%, and the time expected for this event to occur is T ≈ 74 years.These results point out that an optimized slaughter policy is needed, together with investments in technology that improve the reproduction rates and carrying capacity of the Brazilian bovine herd and make possible to reach the population goal in a desirable time interval as T ≈ 12 years with p ≈ 94%, which can be achieved with h 0 ≈ 0.25 and K ≈ 250 million head.

COMPARISON BETWEEN THE MARKOV CHAIN MODEL AND MONTE CARLO SIMULATION
In this study, two indices of viability were calculated: the probability of achieving a particular goal (200 million head) and the mean time for the occurrence of this event.To calculate these indices we used two methods: one based on Monte Carlo simulation, and another using a Markov chain model.To compare these two methods we presented, in Figure 4, both indices calculated with each of the methodologies presented, considering K = 230.Estimates made by the Monte Carlo simulation are in the third column of Table I.An analysis of the curves in Figure 4 shows a great similarity between the two approaches.However, being practical, the modeling by Markov chain requires some careful with the choice of size and the discretization of the space state.These results can be seen in Table III and on the graphs in Figure 4. From the results of Table III we note that, with a maximum slaughter rate of 21%, we reach the population goal with the probability of 63% and a mean time period of 22 years for this event to occur.However, this is not a good policy because the production would be very small.On the other hand, with a maximum slaughter rate larger than 25%, the time intervals necessary to reach the goal become very long, and this slaughter policy is also not advisable with K = 230, but with K = 250 (see Table II).

CONCLUSIONS
In this work, we presented the adjustment of the Richards growth model for the bovine population in Brazil taking into consideration the populational losses due to the slaughter intervention.With the model adjustment that includes the stochastic intervention and posteriorly Monte Carlo simulations, we were able to estimate the indicators of the population viability.By considering deterministic intervention, we calculated the indicators, assuming that the population can be modeled by a Markov chain.We also employed Monte Carlo simulations to compare the results between these two methods.The Markov chain model would be the most natural manner of obtaining estimates of a population viability indicators, but it depends on the discretization of the space of states.On the other hand, the method using Monte Carlo simulation requires the establishment of an upper limit that, in this work, we call T ∞ .This limitation may impair the accurate calculation of the time expected to reach the population goal.
The analysis of the numerical results shows that a slaughter rate higher than 25% can make the growth of the population difficult and, in the future, lead to the reduction of bovine meat production in Brazil.To reach a higher goal, such as 200 million heads in a period less than 20 years, an optimized slaughter policy is fundamental in order to maximize a certain expected return within a finite horizon.
Considering the stochastic intervention, we simulated 10 scenarios changing the maximum slaughter rate and the carrying capacity, and concluded that the goal will be reached in 8 to 20 years (with K = 250).This leads us to believe that the use of a fixed intervention equal to 25% (slaughter in 2005) does not guarantee a sustainable production for K < 240.Given the present scenario, we can conclude that the Brazilian bovine herd will not reach the goal of 200 million heads in the expected time of 15 years.

Fig. 2 -
Fig. 2 -Slaughter of the bovine population in Brazil between 1983 and 2008.

Fig. 3 -
Fig. 3 -Simulation of the trajectories of N t for 30 years ahead.

Fig. 4 -
Fig. 4 -Comparison between the estimates using Markov chain and Monte Carlo simulation.

TABLE I The MLE and SD estimates of the Richards model parameters.
An Acad BrasCienc (2010) 82 (4)