Evaluating Livestock Culling Practices Using Population Projection Models in Neolithic Dalmatia

2.1 Mortality and Fertility Parameters

Fertility and mortality rates were assigned to the set of age classes that can be determined by wear and eruption sequences of caprine mandibular dentition (Payne 1973; Vigne and Helmer 2007; Zeder 2006) and based on a review of relevant literature (see S2). In our simulations, an animal is either slaughtered as part of the idealized or observed culling strategy or dies by other means. Accordingly, culling rates specifically refer to the probability of slaughter, derived from survivorship curves. Theoretical culling strategies are largely based on ethnographic data, but archaeological mortality profiles may include counts of animals that died naturally; Intrinsic mortality therefore specifies susceptibility to disease, starvation, dehydration, predation, genetic disorders, and age-related health deterioration.

The competing influences of culling rates and intrinsic mortality for neonates and infants can have potent effects on herd growth. Increased mortality among lambs and kids is linked to poor environmental conditions (Awemu et al. 1999; Dahl and Hjort 1976; Debele et al. 2011; Debele and Duguma 2013) and some demographic projections assign animals up to 12 months higher mortality rates (e.g., Redding 1981; Upton 1984). However, mortality rates are higher among pre-weaning animals (Hatcher et al. 2008; Mukasa-Mugerwa et al. 2000; Singh et al. 2011; see Table S2.1). The Lefkovitch matrix supports demographic profiles with age classes of unequal length (Crouse et al. 1987; Lefkovitch 1965) so that mortality rates of lambs and kids aged 0–2, 2–6, and 6–12 months can be defined separately. We proportionally allocate lamb and kid mortality into age classes A, B, and C, weighted by the duration of each, with the highest mortality probabilities assigned to the youngest (pre-weaning) animals (Table 1; Table S2.1).

Following Redding (1981), we apply mortality rates of 10% and 15% for adult males and females, respectively. Table 1 gives baseline mortality and fertility rates from which inter-annual variation is simulated in the program (summarized below and detailed in S2). We model culling rates in competition with mortality rates, and assume sick animals die before they are slaughtered (or are slaughtered before they get sick). At each timestep, the survival probability (s) of an individual in age class i is calculated as s_i = 1 – m_i – o_i, where m is the intrinsic mortality rate and o is the culling rate.

Population growth hinges on the introduction of new individuals, which can occur via birth or importation. We exclude importation as a factor contributing to herd growth but acknowledge the importance that animal exchange has in pastoral societies of the past and present. Additionally, the model does not account for delaying the age of first parturition, using contraceptive devices, or manipulating the onset of estrus (Balasse et al. 2023; Rosa and Bryant 2002; Watson and Radford 1960; Wilson 1989). These strategies regulate the timing of the breeding season and affect herd fertility, but, in our simulations, reproduction occurs as a “pulse,” where all births occur at the same point in the annual cycle.

In the simulations herd growth occurs when the number of surviving offspring exceeds the number of animals that were slaughtered or died naturally. Fertility is an estimation of the reproductive potential of a herd, which is dependent on the age of first parturition (set to 1 year, see S2), the number of births per female in a single year (i.e., the parturition rate), and the number of offspring expected per birth (i.e., prolificacy rate). Goats exhibit high twinning rates and shorter gestation periods than sheep, and thus higher prolificacy; peak reproductive performance is attained at 3–4 years for both species (Dahl and Hjort 1976; Erasmus et al. 1985; Mellado et al. 2006; Wilson 1989; Wilson et al. 1984; Wilson and Durkin 1983). The calculation of reproductive parameters is described in detail in S2.

2.2 Survivorship and Culling Strategies

We model the effects of ten idealized caprine culling strategies (Figure 1, Table 2) and five culling profiles from four Neolithic sites on the Dalmatian Coast of Croatia (Triozzi 2024; 2025): Smilčić, Zemunik Donji, Benkovac-Barice, and Graduša Lokve—Islam Grčki (Figure 2, Table 3). Age-structured survivorship probabilities reflect the standardized age class system used by Marom and Bar-Oz (2009).

Figure 1

Figure 2

The zooarchaeological assemblage consists of 238 partial and complete mandibles. Distinguishing between sheep (n = 81) and goat (n = 49) mandibles was aided by morphological criteria (Zeder and Pilaar 2010) and comparative reference materials. Indeterminate caprine mandibles were labeled “ovicaprid” (n = 108; Table S4.1). Based on analysis of wear patterns on 750 teeth, we assigned age groups to each mandible following using a revised version of Payne’s (1973) classification system (Zeder 2006). The revised system distinguishes sheep versus goat wear patterns, and the former schema was applied to mandibles of indeterminate species.

Frequency densities of each age class in the archaeological mortality profiles (Figure 2) were calculated following Brochier (2013). Since tooth wear stages can correspond to multiple age classes, some mandible specimens were assigned to up to four age groups; we proportionally (i.e., uniformly) increased the counts of each possible age group, rather than correcting counts according to the duration of each age class. Survival probabilities for each age class were calculated following Price et al. (2016). To transform survivorship of age classes EF and HI to the standardized age classes we reallocated the probabilities for age class G (4–6 years) and H (6–8 years) using the formula:

where x is the culling rate of each age in age class G (or H), and z is the offtake rate of age class G (or H). This allocation assumes that, for age class G, the same offtake rate applies to individuals aged 4–5 years and individuals 5–6 years. The recalculated survival probabilities are given in Table 3. Survivorship probability p for each age class i is converted to culling rate as o_i = 1 – p_i.

The archaeological mortality profiles we examine are not immune to the interpretive caveats mentioned in the introduction section. Three methodological limitations of the approach must also be emphasized. First, distinguishing sheep from goat skeletal materials based on often subtle morphological landmarks is admittedly less reliable than molecular identification techniques such as ZooMS (Buckley et al. 2010) and can lead to overrepresentation of either species in a faunal assemblage (e.g., Sierra et al. 2023). We acknowledge that a ZooMS analysis would likely shift the proportions of sheep and goats we identified through traditional means. Second, young animals are likely underrepresented in the assemblage. Destruction of immature mandibles by scavengers (e.g., dogs and pigs), trampling, and other taphonomic processes reduce the visibility of caprines and ungulates under 1 year of age (Munson 2000; Munson and Garniewicz 2003). We did not correct for underrepresentation of infants and juveniles (Munson 2000) because pig and dog remains were generally not abundant in the faunal assemblages studied.

Lastly, given the difficulty in assigning sex to caprine mandibles, it is rare for archaeologists to assess past production goals using sex-specific culling rates (cf. Diekmann et al. 2026). Idealized culling strategies typically aggregate male and female slaughtering rates as a matter of convenience and only Payne’s milk, meat, and wool optimization strategies distinguish male from female survivorship, but the rates for females are roughly equal across all three strategies (see Figs 1–3 in Payne 1973). Herders generally view males as expendable and tend to keep herds with higher female to male ratios (Dahl and Hjort 1976; Hadjikoumis 2017; Redding 1981). Assuming males dominate faunal assemblages and the decision to slaughter females is driven primarily by compromised reproductive capacity (e.g., older or barren females; Malher et al. 2001), we apply the idealized and empirical culling rates at face value to males but reduce the female culling rates by 85% (see Section 2.3).

2.3 Population Projections

Herd population dynamics were projected over a 200-year period with initial herd sizes of 150 for goats and sheep. The initial herd size assumes a scale of husbandry similar to Vlach households, where a husband and wife are capable of managing a flock of 300 animals (Nitsiakos 1985: 41), and modern herders on Cyprus, who generally manage herds of 80–150 sheep/goats (Hadjikoumis 2017). Survival probabilities (Figures 1 and 2) simulate the constraints of 15 different sheep and goat culling practices on population growth. An additional “baseline” strategy simulates population dynamics of unmanaged herds. Since herd growth is dependent on females, sparing ewes and does from slaughter would prohibit the comparison of the modeled culling strategies because herd survival would be driven solely by variation in fertility and mortality. Assuming herders value the reproductive capacity of females, we reduced the female culling rate for each strategy by 85%, thereby applying a net female offtake rate of 15%. This adjustment was based on data reporting the culling rate of 14.5% of does in intensively managed herds of goats due to infertility and 17.1% for infertility at 1–2 years (Malher et al. 2001).

Environmental factors affect the quality and quantity of dietary resources and prevalence of disease, which can influence fertility and mortality in sheep and goats (Awemu et al. 1999; Debele et al. 2011; Debele and Duguma 2013; Galina et al. 1995).

At each timestep new sets of fertility and infant mortality rates are sampled from a normal distribution to simulate the effects of environmental variability on herd survival and reproduction (see S2). Sampled prolificacy rates were rearranged so peak prolificacy corresponds to the 3–4 year age class (Erasmus et al. 1985; Wilson 1989; Wilson et al. 1984). At each time-step the mean prolificacy rate (l) was used to calculate the annual reproduction rate (ARR, Equation S2.1). A similar method was used to vary infant and adult mortality rates except that infant mortality decreases with age, and adult mortality increases with age. Simulated variation in fertility and mortality rates is truly stochastic because years with high fertility rates may not necessarily be paired with low mortality rates. The independence of these parameters allows for the influence of herd culling rates on herd demography to be examined with respect to environmental uncertainty.

Demographic transition matrices provide several important metrics including λ (i.e., population growth rate) and the sex ratio of the population. The first metric, λ, is obtained from the dominant eigenvalue of the Lefkovitch transition matrix (see S1) and is synonymous with growth rate: populations increase when λ > 1 and decrease when λ < 1. A steady state is reached when λ = 1. Since fertility and mortality vary at each time-step, λ, sex, and age structure of the population can also vary as a result. Using the initial sets of fertility and mortality rates for 200 timesteps, bootstrapped estimates (n = 1000) of the mean of each λ and the proportion of females were obtained for each herd-strategy combination (Table 4 and Table 5).

2.4 Culling Optimization Algorithm

We applied an optimization procedure to simulate herder responses to environmental impacts on herd growth to explore how flexible culling practices may reduce inter-annual variance in herd growth. The algorithm identifies a factor by which female offtake rates must be altered for the current timestep so that the herd growth rate remains unchanged in the subsequent timestep (i.e., λ = 1). The factor, φ_opt was obtained via an iterative process of minimizing (λ – m)² where m is the multiplication rate of a population in a steady state (Lesnoff 2024). The optimization goal was thus defined as λ = m = 1. Before projecting from t to t + 1, the predicted growth rate, λ, was obtained from the current time-step’s transition matrix. If λ in the current timestep fell outside lower and upper thresholds (i.e., if λ < λ_min or if λ ≥ λ_max) female offtake probabilities were adjusted by φ_opt, the optimized harvest intensity rate. The thresholds λ_min (0.936) and λ_max (1.005) were defined as the 25^th and 75^th percentiles (respectively) of the bootstrapped λ values (λ_boot in Tables 4 and 5), excluding the Baseline strategy. A Levene’s test for equality of variances was performed to determine whether optimizing offtake significantly reduced inter-annual variance in the multiplication rates of sheep and goat herds for each strategy (Table 6).

We elected to use λ rather than herd size to stimulate herder responses for three reasons. The first is a matter of convenience; λ_min and λ_max can be easily changed in the scripting procedure (see S3), allowing for different boundary values to be examined for different culling strategies in a sensitivity analysis (see Section 3.3). Second, we were peripherally interested in whether a stable herd size would emerge in any of the managed populations to understand if some idealized strategies are better suited to larger versus smaller herds. Restricting herd size to an arbitrary number would conflict with our assumption that herders might maximize the number of animals managed by adjusting culling rates. Lastly, herders may be well aware of how many animals are under their care but many management decisions, especially in dairy production systems, are influenced by their perception of the health and fertility of reproductive females (Dahl and Hjort 1976; Halstead 1996; Malher et al. 2001; Talore et al. 2014). This perception is implicitly reflected in the mathematical definition of λ as the dominant eigenvalue of the Lefkovitch population projection matrix. The value λ is paired with a left eigenvector, v, which gives the relative contributions of each stage/age group to long-term herd growth (see S1; Negassa et al. 2015). Therefore, λ subsumes the reproductive capabilities of females in each age class. Using arbitrary thresholds for herd size to initiate a change in culling rates would not account for our assumption that herders understand how herd growth hinges on the reproductive values of different stage/age groups.

3.1 Initial Projections: Herd Structure and Growth

Rapid and exponential population growth is expected for unmanaged herds (Lesnoff et al. 2000; Moritz et al. 2017). The bootstrapped estimated mean λ (λ_boot) for the Baseline strategy is greater than 1 for both goats and sheep (Table 4), confirming our mortality and fertility rates are realistic. If fertility was too low or mortality too high resulting in λ_boot < 1 for the Baseline strategy, populations would decline for all culling strategies. In the Baseline strategy, λ_boot for goats is slightly lower than sheep, contrary to demographic models showing higher intrinsic population growth among goats (Redding 1981; 1984; Upton 1984).

For all theoretical culling strategies, goat herds are predicted to decline (λ_boot < 1, Table 4). The same is true for sheep populations except for Redding’s (1981) Energy and Payne’s (1973) Wool strategies. The lowest λ_boot estimates were attained for the Milk A strategy for both taxa, predicting rapid herd decimation. This is likely explained by very low survival probability of animals under 1 year (Figure 1). Sheep and goat herds are predicted to decline for the Neolithic strategies at Early Neolithic Smilčić, Middle Neolithic Graduša Lokve—Islam Grčki, and Zemunik Donji. Middle Neolithic Smilčić and Benkovac-Barice are predicted to maintain or increase herd sizes (Figure 3).

Figure 3

A significant negative correlation between λ_boot and the proportion of females in the herd (Pearson’s r = –0.847, t = –8.448, df = 28, p < 0.001) suggests high female-to-male ratios are associated with low or negative population growth. This result contrasts with the assumption that female-dominated herds ensure reproductive capacity. Traditional herders keep more females to limit resource consumption by non-reproductive animals (e.g., males and infertile females; Dahl and Hjort 1976; Redding 1981; Wilson 1978). However, the lowest λ_boot is obtained for Milk A, the strategy with the highest proportion of females (Table 4). Low survival probabilities of animals younger than age class of first parturition, set at 12–24 months likely explains the predicted population declines for the Milk A and Early Neolithic Smilčić strategies.

The initial demographic compositions of herds under each strategy are shown in Figure 4. Except for Vigne and Helmer’s (2007) Milk A and Meat A strategies, the number of kids exceeds that of lambs in the first two age classes, but the trend is reversed in subsequent age classes. The opposite of this pattern is shown in the Baseline strategy, revealing the critical role of culling on herd structure. Despite higher prolificacy of breeding does, higher kid mortality leads to fewer goats reaching adulthood. Despite subtle differences in culling rates, this pattern reflects a broader, universal feature of the culling practices examined, which could reflect an implicit tactic of restricting the nutritional demands of herds dominated by adults.

Figure 4

The population projection results are shown in Figure 5. Both taxa decline rapidly under Payne’s (1973) Milk strategy and all of Vigne and Helmer’s (2007) models. Herds decline and stabilize at 50 animals under Redding’s (1981) Security model. Herds grow under Redding’s (1981) Energy and Payne’s (1973) Wool models and are relatively stable for Payne’s (1973) Meat model. For the archaeological strategies only the Middle Neolithic Smilčić and Benkovac-Barice culling rates permit herd growth. The trajectories of sheep and goat populations under the Zemunik Donji and Graduša Lokve—Islam Grčki strategies appear to be stable with smaller herd sizes, and Early Neolithic Smilčić herds are wiped out after 150 years.

Figure 5

To assess whether the initial projections simulated environments that were more hostile for goats than sheep, we replicated the simulations 100 times (Figure 6) and achieved results that were broadly consistent with the original projections. Populations rapidly increase occasionally, perhaps due to concurrently favorable environments but Payne’s (1973) Milk model and Vigne and Helmer’s (2007) models consistently decimate herds when culling rates are held constant from one year to the next. Redding’s (1981) Energy, Payne’s (1973) Meat and Wool, and the Middle Neolithic Smilčić and Benkovac-Barice strategies appear to promote sheep population growth, but goat herd growth tends to be more conservative. Higher goat prolificacy does not appear to compensate for higher kid mortality even in periodically favorable environments.

Figure 6

3.2 Optimizing Culling Rates

Conditionally adapting culling rates prolonged survival of herds (Figure 7). Recall that of all harvest profiles, the relatively low λ_boot estimate for sheep and goats for Milk A predicted the fastest herd decline (Table 4), but this strategy outperformed the other dairy-focused models in terms of herd longevity (Figure 7) and average herd size (Figure 8). Figures 9 and 10 plot 10-year moving averages of herd growth rates for goats and sheep (respectively) with constant and flexible culling rates. The frequency by which culling rates were adjusted is visualized as the degree to which m (green lines) and λ (magenta lines) converge. Convergence indicates similarity between the original culling rates and the optimized ones, interpreted as less intervention. More frequent adjustments of culling rates result in less overlap between the two lines. Additionally, the effect of flexible culling on minimizing inter-annual variance in herd growth rate is observable as reduced distortion of the green (m) relative to magenta (λ) lines. Strategies for which the annual predicted population growth rate (λ) was consistently lower than the threshold λ_min witnessed more frequent adjustment of offtake rates to achieve an actual herd multiplication rate of m = 1. As an example, offtake rates for the Milk A strategy were frequently adjusted to minimize the otherwise rapid decline in herd size. This is visible as relatively smoother green lines that do not converge at all with the magenta lines in the Milk A panels of Figures 9 and 10.

Figure 7

Figure 8

Figure 9

Figure 10

By comparison, herd sizes declined faster for the Meat A and B, Milk B, Fleece (Vigne and Helmer 2007) models and Payne’s (1973) Milk strategy than for the Milk A strategy in the optimization simulation. For these strategies, estimated herd growth rates (λ_boot) were closer to or even greater than the lower threshold (λ_min = 0.936) used to trigger adjustment of offtake rates (Table 4). The optimization procedure was not initiated as frequently because, for these models, λ_min was too low a threshold. As a result, sheep and goat populations were fixed in a perpetual state of decline that was exacerbated by the occasional intensification of culling rates triggered for any time-step where the growth rate exceeded the upper threshold of 1.005. This restriction on herd growth prohibited populations from sufficiently rebounding from previous losses.

The optimization procedure was intended to counteract excessive herd size increases by increasing culling rates when the predicted herd growth rate exceeded the upper threshold. However, this resulted in a long-term negative impact for strategies otherwise amenable to herd growth (i.e., strategies where λ_boot > 1). Despite overall higher growth rates, goat herds declined faster than sheep populations under the Energy, Wool, Benkovac-Barice, and Middle Neolithic Smilčić strategies.

The negative impact on herd size is perhaps most consequential for the Middle Neolithic Smilčić and Benkovac-Barice strategies. Population sizes stabilized during the initial 200-year simulation (Figure 5) and were repeatedly shown to permit modest herd growth under a variety of fertility and mortality conditions (Figure 6). But optimizing offtake caused the gradual decline of both populations (Figure 7). The explanation lies in the suitability of the optimization thresholds used. Predicted growth rates (λ) for these two strategies rarely dropped below the lower threshold (λ_min) of 0.936 (Figures 9 and 10); for most time-steps where λ_min < λ < 1, culling rates were not relaxed, so m = λ = 1. However, more aggressive culling rates were triggered multiple times in the 200-year simulation for these strategies (peaks exceeding the dashed line in Figures 9 and 10) and prevented both goat and sheep herds from rebounding from periods of low fertility and/or high mortality. The sustainability of these strategies inferred from the initial projections is likely attributed to intermittent but nonetheless critical opportunities for population growth that were nullified by the rules of the optimization procedure.

The simulations with flexible culling projected lower average herd sizes for goats than sheep, even for strategies amenable to population growth (e.g., Benkovac-Barice and Middle Neolithic Smilčic; Figure 8). Optimizing sheep offtake rates consistently reduced inter-annual variance in herd growth rates for most strategies (Table 6). In the initial simulation, higher prolificacy and parturition rates enabled goats to outperform sheep in terms of herd size for strategies with growth rates greater than 1. However, in the optimized simulation, if herd growth declines one year but is predicted to increase the next year, culling rates are increased to stagnate growth. The implication points to different long-term impacts of increasing slaughter rates to curb sheep versus goat herd growth. That is, despite overall higher growth rate, inter-annual variability in goat population dynamics cannot be reduced by increasing culling rates alone. This observation also highlights a potential advantage of sheep herding—one afforded by the marginally slower growth rate of sheep than goat populations that made it easier for past herders to minimize variance in the size of flocks of sheep.

3.3 Sensitivity Analysis

We examined how mean herd size and standard deviation is affected by varying upper and lower growth rate triggers for modifying culling rates and female offtake rates in a variance-based (Sobol’) sensitivity analysis. The Sobol’ method (Sobol’ 1990) quantifies (1) first-order indices, which reflect the sensitivity of model results to each parameter in isolation, and (2) total effect indices describing how variance in model results is driven by interaction effects of the parameters. The parameter space consisted of bootstrapped resampling (n = 50) of two sets of 500 values sampled from normal distributions: λ_min ± 0.01, λ_max ± 0.01, and female offtake: 15 ± 5%. Larger parameter spaces are generally preferred but require extensive computing resources, particularly if the analyses covered all 15 strategies for both taxa. Given the redundancies in the results of projections for Vigne and Helmer’s (2007) strategies and similarities between Redding’s (1981) Energy and Payne’s (1973) Wool strategies, we focused the sensitivity analysis only on the three classic, Meat, Milk, and Wool models. The analysis was performed on UBELIX (https://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern.

Figure 11 plots the main effect indices and confidence intervals (95%) for female offtake, λ_min, and λ_max. Female offtake rates responsible for 11.6 to 20.7% of variance in long-term (mean) herd size and 71.6 to 97% of the variance in interannual variations in herd size across all strategies (Table 7). Figure 12 illustrates, perhaps unsurprisingly, the dominance of this parameter’s interaction with λ_min, and λ_max. These results suggest that herd longevity and stable production are vulnerable primarily to the survival of females, a finding that generally conforms to ethnographic studies of traditional pastoral systems (Dahl and Hjort 1976; Dyson-Hudson and Dyson-Hudson 1969; Galvin 1985; Hadjikoumis 2017; Halstead 1998b; Nitsiakos 1985). Adapting culling rates only marginally impacts herd survival when reproductive animals are not managed carefully.

Figure 11

Figure 12