A comprehensive transition matrix model for projecting production and resource consumption in reindeer herds

A deterministic herd model was developed for projecting the dynamic changes in reindeer herd size and structure under defined harvest policies. The model distinguishes between females, males and castrates up to an optional number of age-classes. Calves are further classified based on the age and status (present/absent) of their mother. The yearly cycle is divided up into a maximum of 11 time steps, including five grazing seasons. The model is described in general terms using the Leslie matrix approach in order to suit different computer implementations. The conventional Leslie matrix solution was extended so that nonlinear features and stochastic variation in performance parameters could be considered. Computational procedures for making detailed economic evaluations of harvest output and herd feed requirements or consumption are given. This general purpose model can be tailored to specific study conditions. A n advantage of this is that the sensitivity to necessary approximations can be tested with the general purpose model. The model is intended for use in both research and extension work. K e y w o r d s : populat ion dynamics, herd simulat ion model , compartment model , Leslie matrix approach, deterministic model , event-based model , stochastic environment, nonlinear features, populat ion control , harvest evaluation, product ion resource consumption. Rangifer, 14 (3): 99-112


Introduction
Reindeer husbandry is practiced in northern Scandinavia and Finland, with the main aim being meat production.The economic return is determined by the amount as well as the quality of meat produced.Since production is largely based on naturally produced forage in the boreal forest and mountain areas, the pasture resources must be sustainably managed.There is little reliance on supplemental feeding except in the southern part of the reindeer husbandry area in Finland and the parts of Scandinavia, where the pastures were contaminated with radioactive fallout from the Chernobyl accident.A reindeer production system can only be efficiently managed if it is possible to foresee the demo-graphic effects of particular harvest policies and management-induced changes in the external conditions.Like most biological production systems this system involves interactions between system components as well as various time lags and feedback mechanisms.As result of this great complexity, reliable predictions are very difficult to make.
Simulations based on dynamic models can contribute substantially to our understanding of systems behaviour in wildlife management and ecological research.In reindeer production this approach has not been used much till now, probably because knowledge about biological production parameters, such as animal performance under different conditions, pasture productivity, etc., has been very limited.Extensive amounts of production data, for use in model calculations, have been collected in Norway since the 1970s (Lenvik, 1988) and in one Saami village herd in Sweden since 1985 (Petersson et al, 1990).The use of dynamic modelling as a tool in reindeer research has been discussed by Danell (1985) and Eberhardt (1991), among others.Petersson & Danell (1992) described an application of an earlier attempt to model reindeer production systems and recent attempts were described by Virtala (1992) and Moxnes et al. (1993), all using different techniques.The firstmentioned study focussed only on the reindeer herd, while in the other two both animals and pasture (lichen) were included in the model.Model calculations, although not formalised with a unified model, were reported by Lenvik (1988).
The aim of this paper is to present a simulation model for use in predicting reindeer herd productivity and resource consumption.The model is comprehensive in the sense that there are linkages between population and management submodels.The intention was to develop a general computational tool for use in solving a wide range of management and extension problems in reindeer husbandry.This model differs from most others developed for use in wildlife management and ecology in that a more detailed information about population demography and animal performance is presumed than usually is possible for wild animals.The problem of identifying relevant paramenters for the model was dealt with separately (Petersson & Danell, 6., 1993).

System delineation
The purpose of this section is to describe the reindeer production system at a level of resolution that is commensurate with our modelling objectives.On the basis of this description, the approximations, which are necessary in the development of the dynamic herd model, become more definable.Key characteristics of the production system are summarized in Fig. 1.For modelling purpose it seems most convienent to consider the reindeer herd and the pasture resources as separate, but interacting, subsystems owing to their completely separate biological frameworks and different time scales with regard to dynamic behaviour.Here, the herd is viewed as the nucleus of this conceptual system.Pasture resources are, however, an equally impor-tant part of the current production system.Clearly, the long-term abundance of natural foods determines the overall limits on production, whereas variation in the abundance and availability of forage will influence the dynamic structure and productivity of the herd also in the short term.In view of this interaction, long-term optimization of the reindeer production system necessitates the simultaneous recognition of both subsystems.When making short-term projections it might be appropriate to treat the availability of pasture resources as an external factor.
The dynamics and productivity of the herd subsystem are influenced by a number of additional factors, as is indicated in the upper right-hand parts of Fig. 1.Management measures might, however, be corrections directly caused by the development within the herd subsystem or indirectly via the effects of the animal herd on the pasture subsystem.In this context we consider all management measures, except those steering the harvest, as external conditions for the animal subsystem.
The economic return depends on revenue based on the quantity and quality of meat produced and its market value, and the production costs.The latter consists partly of base expenditures, connected with the type of production technique, and partly of e.g.expenditures for supplementary feeds or the treatment of individual animals, which are determined by the structure and development of the herd subsystem.Other measures of system properties might also be useful, e.g.sensitivity of production to external influences or long-term efficiency with respect to the utilization of pasture resources.
The herd subsystem is perceived as a dynamic system, in which calves are born each year into the subsystem, and animals within each sex cohort become successively older and ultimately pass out of the system as a result of death losses (e.g.accidents, predation, disease, starvation) or harvest.Castrates are an optional intermediate stage in the transformation of males to slaughter animals.Calf crops are determined by the number and reproductive performances of the females.Females also influence the survival and growth of their calves.In general, survival rates, reproductive and mothering abilities, and individual weight changes are related to sex, age and external factors, and are also affected by stochastic variation.The adult sex ratio and age structure of males may also influence the reproductive rates of the various female age classes (Lenvik et al., 1988).Likewise the quality (or commercial value) of carcasses and consumption of pasture and other production means are influenced by the sex and age structure of the herd.

The production year
Detailed descriptions of the annual cycle in reindeer production systems can be found elsewhere (e.g.Skjenneberg & Slagsvold, 1968;Skjenneberg, 1989).There are obvious seasonal differences during the production year as concerns animal physiology and performance, the type and avaliability of feeds, and other external environmental conditions, management activities, etc.For modelling purposes it seems more advantageous to divide the production year up into discrete and homogeneous periods, for which typical conditions can be defined, than to consider it as a continuous period of time.Alternatively one can divide the year in periods delineated by 'natural' break dates.Such dates could be points of time at which animal numbers (births, slaughters) change abruptly or management activities enable animal numbers and performance to be recorded.
Relevant major events during the production year under conditions typical for reindeer husbandry in Fennoscandia are summarized in Fig. 2. Generally, four major grazing periods can be distinguished, which differ in location of the grazing area, forage quality and environmental conditions in general.None of these periods are homogeneous in terms of environmental conditions, nor can they be distinctly separated.For example, early and late halves of the 'autumn', 'winter' and 'spring' grazing seasons, especially, may differ considerably in forage type and availability as well as weather conditions.
In our situation, it was most convenient to consider calving, corralling at calf marking and the two slaughter periods as major breaks dividing the year into four major 'grazing periods'.However, these breaks are actually periods rather than distinct dates.Due to the often dramatic variations in winter conditions and the long period from autumn slaughter to calving, a further subdivision of the 'winter' period can be justified.Consequently, up to five discrete periods, in addition to birth and slaughter periods can be a proper level of detail in discribing the yearly cycle.or key events as described in Fig. 2

Basic outline
The basic approach is to approximate the structure and dynamics of a reindeer herd using a compartment model with separate compartments for calves (age < 1 year) of each sex and dam age as well as for each age-class of adult (age > 1 year) females, males and castrates.Motherless calves are assigned to separate compartments according to their sex and the period of the year during which the dam was lost.
Time is advanced in, at most, 11 discrete steps per year, corresponding to the sequence of periods and events in Fig. 3.Each time step is associated with specific processes as desribed in Table 1.The six pulses representing aging, births and harvests, are assumed to occur at the midpoints of the periods they mimic.Therefore, periods Pi -P5 are presumed to include adjacent parts of the periods modelled by pulses.For biological reasons, aging and births may be considered as reasonably fixed in time.Other 'break dates' are varying in time due to the management or specific circumstances of the study situation (different production techniques, regions, etc.).

Transition equations
Basic dynamic processes include the generation of newborn animals into the modelled system via births, the successive transition between compartments at aging, castration of males, loss of maternity, and flow out from the modelled system via nonslaughter-related death losses and harvest.Regular animal types nxi for x = 1 (female calves), 2 (male calves), 3 (females), and i = 1, ...k3 (= max.considered female age in years) nxi for x = 4 (males), 5 (castrates), and i = 1, ...k4 or ks (= max.age in years for each type) Motherless calves mxp and x = 1 (male calves), 2 (female calves), for p = 1, (mother lost in Pi), ... 5 (ditto in P5) The model is event-based with time steps that tional convenience are organised like a Leslie mathave not even width.It is formally defined as a set rix (Leslie, 1945;1948).The Leslie matrix format of linear difference equations, which for computa-has commonly been used in studies of population Rangifer, 14 (3), 1994 dynamics and a wide range of other linear transition problems by other authors.This does not, however, preclude its application to nonlinear behaviour since model parameters can also be modified during computations according the state in certain compartments, e.g.density dependent survival probabilities or calving probabilities depending on the adult sex ratios during the rut.
Let nxi(t-l) and mxp(t-l) be state variables (abundance of animals) for the various compartments, as defined in Table 2, evaluated at the end of the previous time step (f-1).The passages of animals over time can then be expressed by simple linear transition equations as follows: During time steps corresponding to periods Pi to P5, only surviving animals are transferred; these appear in the state variable at the end of time step t.For adult animals and already motherless calves the transition equations are and where sxi and are the survival probabilities of animals in subclasses xi for x = 3, 4 and 5, and xp for x -\ and 2, respectively, during this particular time step.When females with a calf are lost, similar proportions (i.e. 1 -s3l) of surviving calves in the regular calf compartments of each sex and corresponding to the particular age classes of lost mothers are transferred to the appropriate compartments of motherless calves according to sex of calves (x) and the period (p) when the mother was lost.Thus for the regular calf compartments (x = 1 and 2) we have nx,(t) = s3l sxl nxl{t-l) (3] and for the compartment of calves becoming motherless in the present time step State variable mxp is zero, as a result of the aging in time step A (see below), until that the transfer given in [4] has occurred.
Survival probabilities operate on the animal numbers at the beginning of the time step.Therefore they are dependent on the step duration.Since period lengths are not explicitly fixed, the probabilities need to be tailored to the actual period lengths.This can be done via the intrinsic or daily survival rate, rxl, as where d is the actual period length in days, assuming a constant rxi over the period.Appropriate rate values can be obtained by combining equations (1) and [5], given that the abundances (e.g.nxt) at beginning and end of a reference period and the period length (ct) are known: In the slaughter steps, the retained animals are transferred from animal states at time t -1 to the states at time t in the same manner as in (1), i.e. nxl{t) = gxlnxl{t-\), ( 7) in the regular animal compartments, and for motherless calves, where gxi, g^ and hxt, hxp are the proportions retained at the involuntary and intentional slaughter steps, respectively.Involuntary slaughter (different from natural mortality where carcasses are not utilised) is defined as culling due to an animal's poor condition, an injury, etc., whereas intentional slaughter refers to harvesting according to a chosen harvest policy.The latter is the measure by which herd size, age structure, product output, etc., is regulated.Thus, hxi and hxp are functions of animal abundance and require the evaluation of nxi (f-1) and mxp (t-1) prior to performing ( 8) and (10).Management manipulation of slaughter for regulation of standing stock to an intended size or structure introduces a nonlinear feature to the model equations.
The proportions retained at intentional slaughter are in principle calculated as where h x l (t) is the intended number of individual to retain in age class xi.This number can be derived in various ad hoc ways.The way chosen here for steering the slaughter is to describe the intended slaughter policy by defining the maximum herd size or number of animals in each sex or animal type to be retained, a priority order for harvest and the minimum proportion to be retained in each animal age class.The highest priority animals are harvested first, followed by the next highest until the inteded maximum number of animals is reached.Classes with the same priority are reduced in parallel till the maximum scope for slaughter is reached in these classes.
Male castrations, when done, are usually carried out in connection with summer corralling and prerut slaughter.In the model this is equivalent to transferring defined portions (q ) of males to the compartments of castrates of the same age, i.e. nSl(t) = q n4l (t-1) (12) and reducing the male states accordingly.
At the aging time step (A) animals are passed with a unit probability to either the subsequent age class, or, if i represents the last considered age class, to the same age class a time step later.For adults (x > 3) this is and where &x denotes the last age class of animal type x.
Aging of calves is performed by collecting them into the lowest age classes of females and males (x = 3 and 4), as and equating all «x, (t) and mxp (t) for x = 1 and 2 to zero. In where k3, k4 and ks are the highest considered age classes of females, males and castrates, respectively.
A general expression for the transitions performed in one time step is then Rangifer, 14 (3), 1994 where P(î) is an extended version of the transition probability matrix defined by Leslie (1945).P is a square matrix, partitioned to conform with the structure of n: The 25 blocks of P are involved in the passage of animals as indicated in Fig. 4. A different P is required for each of the 11 time steps during a simulated production year.
For time steps Pi, P2, P3, P4 and P5, P is blockdiagonal, i.e. non-zero values occur only in the diagonal blocks and all off-diagonal blocks are null matrices.The blocks P33, P44 and P55 are diagonal matrices with survival probabilities of the different age classes of adult animal types for the particular time step (not labeled specifically here) as elements, e.g. for females calves.The diagonal elements of the lower diagonal block execute similarily equation (3).Likewise, the non-zero elements in the upper off-diagonal block of Pn sum the surviving motherless calves from each dam-age compartment of calves (corresponding to the column numbers in Pn) into the appropriate compartment of motherless calves (corresponding to the row number in Pn), as defined in (4).
In the PSL1, PSL2, ASL1 and ASL2 steps, P is a diagonal matrix with the proportions retained at slaughter, as defined in ( 7) to (10), in the diagonal positions.
The transfer of males to castrate compartments is acomplished by entering the transfer coefficients (q) from ( 12) on diagonal entries of P54, 1 -q on diagonal entries of P44 and entering zero elsewhere.This can be done in a separate iteration of ( 18) with all other blocks of P kept as null matrices or, simultaneously, with the other transfers in, for example, the P4, PSL1, P5 or ASL1 steps by modifying P for these steps accordingly.
The aging of adults in time step A is accomplished by putting unit probabilities of ( 13) and ( 14) on the subdiagonal and the last diagonal positions

/
The corresponding elements are entered in P23 to achieve male births.
The iterations of ( 18) can be done either in individual time steps or in blocks of time steps which do not require separate evaluations of n.Occasions when n must be explicitly known are determinations of intentional slaughter for regulation of herd size and evaluations of harvested animals.A minimum set of computations for an annual cycle is shown in Table 3.However, the modification of the fertility coefficients in the transition probability matrix of the birth step based on the adult sex and age structure at rut and the computation of pasture and feed consumption and production costs based on herd structure would require the first computational step in the table to be split in several separate steps.

Projection of harvest and slaughter revenues
As shown in Fig. 1, meat and the revenues thereof are considered as an output from the herd system.
The number of harvested animals is found by summation of (! ~ gxi) nxi (£-1), (1 -gxp ) mxp{t-l) and (1 -hi ) na{t-\), (1 -hxp) m xp (t~l) in the two adjacent slaughter steps at pre-rut and autumn slaughters, respectively.For one slaughter step (from £-1 to £), these summations are in matrix notation: compartmentwise or in monetary units Rangifer, 14 (3), 1994 where a; is a possible subsidy price per slaughtered animal, /3; is a slaughter and processing cost per animal, and i)j is an average carcass vahie computed as a weighted mean of animals valued differently per kg weight (possible per kg subsidy plus market pri-ce) in different carcass weight classes for; = 1, size of n.
For non-selective harvest within age groups, can be found as follows, assuming that weights are normally distributed in each animal type-age group.Let wjftj be a mean weight (elementof w'^) and Oj(t) be the standard deviation of individual weights of the particular animal type and age group.Furthermore, let T0, r1; r2,... T^, ... be the fixed weight class borders, as exemplified in Fig. 5  where $>( n ~wm ) ;s tne standard normal distribu¬ tion function giving the proportion between -co and value of class border k. and <j > ( n ~ ^ ) is the standard normal density function evaluated at class border k.The former is an integral which has to be approximated numerically with standard routines, and the latter is explicitly given by the wellknown expression for the normal distribution density.Now, the number of animals from compartment j found in weight class k is simply p^ times the harvested number in this compartment (element ;' of n./, (t)), and the mean weight of the animals found in this class becomes

+ &jk a j(t)-
The average carcass value #; is then the weighted mean of times the price per kg for the different price classes, using the corresponding pjk as weights.

Pasture and feed consumption
The herd model can be connected to pasture and feed supply models via information on feed requirements or feed consumption computed from animal abundance in different periods.The simplest alternative is to define a vector y^ which includes the average requirement or consumption (y^t)) of pasture or other feeds per animal in type-age class ;' during each period of interest (f).As feed requirement or consumption usually is not directly proportional to animal weight, can be found as a weighted mean based on a function f(w) relating consumption to animal weight (assumed normally distributed), as y]{t) = \l\ <rM fW dw ^ j Z\ fH dw where <( > (w) is the normal density function of the animal weights in type-age class; at time t, with w --N(w;^, ofa ) and wx and w2 chosen symmetrically around w'^ty A possible form of f(w) is for example where d is period duration in days and z is the daily energy requirement per kg metabolic body weight, w 0 J i , as applied by Lenvik (1985).The integral expressions can easily be solved by numerical integration.The total consumption during period t is then [32)

Including stochastic effects
Although the modelling approach is deterministic, stochasticity in the paramenters can be included in the P matrices to account for stochastic environments.Individual animal variation can also be introduced by drawing the parameters from the appropriate distributions rather than using deterministic mean values.
Thus, stochastic survival probabilities, sxt^, can be found as ( r ) where X is a random deviate, common for all animal classes (disregarding the nonlinearity of probability values), drawn from a distribution typical for average survival probabilities sxi^ over time (e.g.different years) accounting for stochasticity in average survival level, and KXI is a random binomial deviate generated numerically from a binomial distribution for samples of size nxl(t-\) and mean probability s xi(t) + ^ (eg- Press et al., 1989 p. 208-209) accounting for the sampling noise in each animal class.Stochastic fertility coefficients can be generated in a similar way.Thus, both require evaluation of the actual abundances at the end of the previous time step.
Stochastic average weights, wjfy, for a compartment with nxi animals can be computed where i\x (common for all animal classes) and r\2xl are randomly generated standard normal deviates (e.g.Press et al, 1989 p. 204), cr5>, is the standard deviation of mean weights between years and am, is the within-year standard deviation of the individual weights for the particular animal type and age, assuming normal distributions of mean weights between years and individual weights within year for each animal type and age.Thus, the first stochastic term accounts for the variation between years and the second for the sampling variation as depending of sample size.Other distributions than the normal are also possible here but more difficult to handle computationally.
In all of the above cases, one or both stochastic components could be included.With none included, (33) and (34) reduce to the original deterministic values.

Introduction of nonlinear properties
Besides harvest for regulation of herd size, a few other nonlinear features may easily be included in the model computations.Examples are densitydependent survivals, fertilities or animal weights, and fertilities dependent on male frequency and age distribution.These involve modification of the de-terministic probability and animal weight parameters in relation to the total animal abundance viz-aviz the available feed resources and adult male to female ratios, respectively.Different ad hoc solutions are possible, provided the appropriate causal relationships can be assumed.
It is also possible to consider selection effects on animal weights, for example in (26), if the selection applied can be defined.For truncation selection (all animals above or below a certain weight), standard expressions for computing mean weights of the slaughtered as well as the retained animals are available (e.g.Falconer, 1989).When, for example, the lightest animals are harvested, the elements of W( £ ) are where D is the standardised selection differential found as -<K^-)/(lhj [t]) (e.g.Falconer, 1989) for the truncation weight T and the proporti-

Discussion
The model presented in this paper is an extended version of the one used by Petersson & Danell (1992) for estimating herd production losses due accidental losses of animals.This earlier model was developed in mid-1980s and used in extension work to motivate recording of performance data.The principal structure was similar, but it was formulated as a 'flow chart' suitable for coding in the DYNAMO simulation language (e.g.Gustafsson, 1983).
The structure of both models differ in complexity from the models found in the literature for studies of population management and ecology.Usually the sex and age structure is simplified so as to only include a few age compartments.Most discrete time models use one reproduction cycle or one year as time step.Furthermore, few population models have been developed for projection and analysis of reindeer herd dynamics.Two recent examples are models, that were developed for studies of the reindeer-lichen system with focus on the current overgrazing of winter pastures in the Fennoscandia area.The herd model used by Virtala (1991) is a discrete time model with no age structure but with four periods per year as time steps ('winter', 'calving', 'summer' and 'harvest'). Moxnes etal. (1993) used a time-continuous model with three age classes in each sex held separated (calves, yearlings and older).The reason for choosing a time-continuous model was to facilitate formal representation of nonlinear relationships.In our model this is largely solved by splitting the year into several time steps, which allow a detailed monitoring of dynamic changes in animal numbers within years.The consequences thereof can be considered during in the following time steps by changing the probabilities in the transition matrices accordingly.Interactions with pasture abundance can be treated similarly.
In this paper we suggested a fairly comprehensive model for simulating reindeer herd dynamics.The purpose is twofold: The first objective is to provide a common basis for studying a wide range of dynamic problems in reindeer husbandry.The advantage of this, we think, is a gain in comparability by being able to apply a unified model to different studies.We also pointed out various possibilities for further refinements, which consider non-linearity and variable external conditions.This was done mainly to ensure that the approach does not become to limited in its application.
A second reason for outlining a very extensive model is to elucidate the extent and level of detail of knowledge, that might be required in various study situations.In particular situations, the degree of approximation can be quantified by testing the model for sensivity to errors in incompletely known input parameters.During the development of the model we considered assessing its sensitivity to variation in parameter values in general terms.However, we concluded that such work would have little meaning because of possible parameter interactions and the difficulties in defining a result function, which is generally adequate.
An important disadvantage with a complicated model is that irrelevant or minor factors included in the model might contribute to noise, and thereby render it more difficult to focus on major factors and arrive at general conclusions about them.We therefore recommend that one should avoid differentiating parameters too much with regard to age or season if solid data are lacking.An alternative would be to test the model's sensitivity to detailed differentiation in the particular situation studied.
Rangiier, 14 (3), 1994 Even after the model has been coded into a computer software programme, its operation can still be simplified when adequate input data are lacking, by aggregating or simplifying the model parameterization.For example, one can use unity survival probabilities for time steps P2 and P3 to mimic a yearly cycle with only three grazing periods.If one slaughter occasion (e.g. the PSL1 and PSL2 steps) is dropped by steering the model to zero slaughter that period, the autumn time step P5 can also be merged into the adjacent period in order to get only 'winter', 'calving', 'summer' and 'harvest' steps in the model.Likewise it is possible to ignore the differential survival rates of calves with dams of different age or mother status (i.e. has mother vs motherless), use the same performance parameters for all higher age classes of males and females and set castration probabilities to zero in order to get a two-sex model with only a few age-classes in each sex.
When applying the model there will obviously be problems in finding consistent parameters.A set of conceptual 'sub-models', which compile the knowledge about animal performance levels and variation between seasons of the year, different years, regions, production types, etc., would be helpful in preparing appropriate input parameters.This demands considerable research within the area of reindeer biology and husbandry.Work in this area has been reviewed and reported elsewhere (Petersson & Danell, O, 1993;Petersson & Danell, B., 1993 a, b).Records on production parameters have been collected in herds with individually tagged animals in Norway since the 1970s (Lenvik, 1988) and in Sweden since the mid-1980s (Danell, 1986;Petersson et al. 1990).
The transition matrix approach is by far the most commonly used technique for population modelling, however there are many variants on the general approach.It has the advantage of having a welldefinied structure, and computations are easily understood since they are based on standard matrix operations.Matrix representation of compartment models gives access to a considerable number of mathematical tools for model analysis and optimisation of population dynamic problems (see e.g.Cullen (1985), Caswell (1989) and Getz & Haight (1989) for reviews).The complexity of the suggested model prohibits the use of most analytic techniques, however.When analytical methods are required for solving particular problems, a useful procedure might be to first simplify the model considerably and try to derive general (and maybe qualitative) results to serve in guiding further studies with a more complete model.The method when using the latter approach would be to survey interesting parameter combinations by simulation.Prior knowledge of general results would most likely help to limit the range of parameter combinations requiring screening and simplify the search for critical points, which is difficult when many parameters are involved.With too many parameters one is forced to deal with a huge number of parameter value combinations or to examine fewer levels of each parameter.
The described modelling approach fits a number of different implementation tools.Examples include specialized simulation languages such as DYNAMO (e.g.Gustafsson et al, 1983;Petersson & Danell, 1992) or DEMOS (e.g.Gustafsson & Liss, 1985), matrix software such as the interactive MATLAB (MathWorks, 1989), the SAS interactive packages for matrix computations (SAS Institute Inc., 1990) or commonly available spreadsheet programmes.Our choice has been FORTRAN, owing to the flexibility allowed in the use of shortcuts when handling large sparce matrices, the ease in tailoring the program to specific study situations or the use of specialized numerical routines, and the portability of programs to different computer hardware.Considerable computing efforts can be saved, for example, by treating diagonal matrices as vectors and modifying the matrix computations accordingly.
From computational point of view, the use of a large, complicated model can be disadvantageous.In the past, simulation models were often kept simple to reduce the computational burden.However, this is no longer a concern today, even for simulations covering periods of several decades.

Fig. 1 .
Fig. 1.Description of the reindeer production system, with emphasis on the influences of reindeer herd dynamics on productivity.The subsystem is delimited by the dashed lines.Solid arrows represent direct effects of one component on another, whereas dashed arrows indicate possible influences, which are considered indirect in this context.Broken bold arrows within the herd box represent transitions in the status of animals over time.

Fig. 2 .
Fig. 2. Main characteristics of the reindeer production year.

Fig. 3 .
Fig. 3. Simulation time steps corresponding to periodsor key events as described in Fig.2and Table1.Broken arrows indicate 'break dates' which may vary according to external conditions.
the subsequent time step (B) new animals are generated into the regular calf compartments via births determined by fertility parameters defined for each female age class.Let /J be the probability that a female of age i has a calf which survives the neonatal period, and ux the proportion of each sex among born calves.Then we have n xi (*) = u xf,n3i{t-l)for x = 1 and 2. (16)Matrix representation of transition equationsThe above transition equations exceed 500 per year cycle if 10 or more age classes of adults are considered.In computer implementation it is convienent to handle the computations in matrix terms.Let n' be the transpose of column vector n with elements nxi and mxp organized as

Fig. 4 .
Fig. 4. Involvement of blocks of the transition probability matrix P in different transition processes.
P 2 2 become more complicated because motherless calves are involved.For example, for female calves in time step P4 we have (equation male calves in P22.(The survival probabilities in parentheses in the upper left block are not needed in this particular time step since the corresponding elements of n (t-l) are still null.)When performing the multiplication in (18), the diagonal elements of the upper diagonal block execute equation (2) for all compartments of transferred to the first age class of adult females and males by letting the rectangular P31 and P42 have unit transfer probabilities on the first row and zeros elsewhere.All other blocks of P are null matrices in this step.Births are accomplished in step B by introducing the transfer coefficients for births from (16) on the diagonal proceeding from the lower right-hand corner of P]3 and P23 and again leaving all other blocks of P as null matrices.Thus, to achieve female births A minimum set of computations required to fulfill a yearly cycle including projection of slaughter revenues.Matrix computations are iterations of equation(18]  given in the text.Subscripts of P refer to time steps as defined in Fig.3, is the product of 6 transition probability matrices.n(end P4) = P(P4) P(P3) P(B) P(A) P(P2) P(P1) n(end ASL2) = P(I) n(end ASL2) Projection of PSL1 harvest from n(end P4) n(end PSL1) = P(PSL1) n(end P4) Evaluation of n(end PSL1) to determine PSL2 slaughter Projection of PSL2 harvest from n(end PSLl) n(end PSL2) = P(PSL2) n(end PSLl) n(end P5) = P(P5) n(end PSL2) Projection of ASL1 harvest from n(end P5) n(end ASL1) = P(ASL1) n(end P5) Evaluation of n(end ASL1) to determine ASL2 slaughter Project of ASL2 harvest from n(end ASL1) n(end ASL2) = P(ASL2) n(end ASL1)

Fig. 5 .
Fig. 5.The distribution of individual weights of animals in compartment j at time t divided in three weight classes representing different carcass or meat values per kg weight.Symbols are explained in the text.108 on 1 -hfo) harvested (element of H(t )).This assumes that weights are normally distributed, and the selection criteria are the weights themselves.The mean weight of the retained animals are found by substituting § (-^2*»-) / hm for D in (35).With selective harvesting, the distributions of weights in the harvested and retained groups are nonnormal, however, which interferes with the normality assumptions necessary for (29) -(31).

Table 1 .
Description of the 11 simulation time steps deviding ; up the year (Secondary processes in parentheses)..

Table 2 .
Definition of state variables displaying the abundance of animals according to animal type and age.