Table 1
Parameters used to vary fertility and intrinsic mortality for sheep and goat populations. Prolificacy (litter size) is based on mean for sheep (1.15 ± 0.122) and goats (1.49 ± 0.275). Parturition rates assume a 300-day inter-birth interval (see S2). Adult mortality is based on Redding (1981). Infant mortality assumes greater susceptibility of pre-weaning animals to disease and starvation. Parturition and prolificacy rates shown here are stochastic; the model varies mortality based on rates above.
| AGE (YEARS) | PARTURITION | PROLIFICACY | MORTALITY | |
|---|---|---|---|---|
| FEMALE | MALE | |||
| Taxon: Goat | ||||
| 0.17 | – | – | 0.225 | 0.225 |
| 0.50 | – | – | 0.150 | 0.150 |
| 1.00 | – | – | 0.075 | 0.075 |
| 2.00 | 1.23 | 1.01 | 0.125 | 0.150 |
| 3.00 | 1.80 | 1.48 | 0.125 | 0.050 |
| 4.00 | 2.42 | 1.99 | 0.075 | 0.050 |
| 6.00 | 2.46 | 2.02 | 0.075 | 0.050 |
| 8.00 | 2.05 | 1.69 | 0.150 | 0.500 |
| 10.00 | 1.57 | 1.29 | 0.500 | 0.500 |
| Taxon: Sheep | ||||
| 0.17 | – | – | 0.160 | 0.160 |
| 0.50 | – | – | 0.107 | 0.107 |
| 1.00 | – | – | 0.053 | 0.053 |
| 2.00 | 1.24 | 1.02 | 0.125 | 0.100 |
| 3.00 | 1.45 | 1.19 | 0.125 | 0.050 |
| 4.00 | 1.55 | 1.27 | 0.075 | 0.050 |
| 6.00 | 1.60 | 1.32 | 0.075 | 0.050 |
| 8.00 | 1.51 | 1.24 | 0.150 | 0.500 |
| 10.00 | 1.43 | 1.17 | 0.500 | 0.500 |

Figure 1
Survivorship curves from ten theoretical culling strategies.

Figure 2
Mortality profiles for four Neolithic sites on the Dalmatian Coast of Croatia (Triozzi 2024), with calculated survivorship curves.
Table 2
Percentage of herd predicted to survive under ten theoretical culling strategies associated with different production strategies standardized by Marom and Bar-Oz (2009). Survival probabilities derived from aRedding (1981), bPayne (1973), and cVigne and Helmer (2007).
| AGE CLASS | ENERGYa | SECURITYa | MEATb | MILKb | WOOLb | MEAT Ac | MEAT Bc | MILK Ac | MILK Bc | FLEECEc |
|---|---|---|---|---|---|---|---|---|---|---|
| A (0–2m) | 90.4 | 90.4 | 85 | 47 | 85 | 81 | 86 | 22 | 83 | 69 |
| B (2–6m) | 90.4 | 90.4 | 75 | 42 | 75 | 34 | 68 | 11 | 50 | 35 |
| C (6–12m) | 77.6 | 64.5 | 70 | 39 | 65 | 11 | 28 | 4 | 36 | 24 |
| D (1–2y) | 47.6 | 38.0 | 50 | 35 | 63 | 7 | 6 | 3 | 18 | 17 |
| E (2–3y) | 25.0 | 25.0 | 30 | 28 | 57 | 7 | 6 | 3 | 18 | 17 |
| F (3–4y) | 23.9 | 23.9 | 22 | 23 | 50 | 3 | 1 | 2 | 6 | 6 |
| G (4–6y) | 18.2 | 18.2 | 19 | 18 | 43 | 1 | 1 | 1 | 1 | 1 |
| H (6–8y) | 16.1 | 16.1 | 19 | 18 | 43 | 1 | 1 | 1 | 1 | 1 |
| I (>8y) | 11.8 | 11.8 | 10 | 10 | 20 | 1 | 1 | 1 | 1 | 1 |
Table 3
Survivorship probabilities derived from age-at-death data associated with four Neolithic sites examined. EN = Early Neolithic; MN = Middle Neolithic. Probabilities calculated following Price et al. (2016).
| AGE CLASS | BENKOVAC-BARICE MN | GRADUŠA LOKVE—ISLAM GRČKI MN | SMILČIĆ EN | SMILČIĆ MN | ZEMUNIK DONJI MN |
|---|---|---|---|---|---|
| A (0–2m) | 98.61 | 81.25 | 100.00 | 98.25 | 100.00 |
| B (2–6m) | 88.43 | 81.25 | 72.73 | 90.64 | 82.83 |
| C (6–12m) | 85.88 | 67.19 | 56.82 | 83.92 | 70.20 |
| D (1–2y) | 54.86 | 37.50 | 19.32 | 41.23 | 27.27 |
| E (2–3y) | 39.58 | 28.12 | 13.64 | 37.43 | 21.21 |
| F (3–4y) | 28.70 | 21.88 | 6.82 | 27.92 | 3.03 |
| G (4–6y) | 1.63 | 4.80 | 1.14 | 4.95 | 0.00 |
| H (6–8y) | 1.63 | 4.80 | 1.14 | 4.95 | 0.00 |
| I (>8y) | 0.00 | 3.18 | 0.00 | 1.10 | 0.00 |
Table 4
Bootstrapped (n = 1000) means of predicted herd growth rate, λ (λboot), proportions of female and male sheep and goat in herds when survivorship probabilities derived from theoretical culling profiles (Payne 1973; Redding 1981; Vigne and Helmer 2007; Table 2) are parameterized as offtake rates in the Lefkovitch population projection matrix (Lefkovitch 1965). Values shown here reflect 200 years of variation in mortality and fertility rates given in Table 1.
| ENERGY | SECURITY | MEAT | MILK | WOOL | MEAT A | MEAT B | MILK A | MILK B | FLEECE | BASELINE | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Taxon: goat | |||||||||||
| Lambda | 0.992 | 0.985 | 0.989 | 0.950 | 0.996 | 0.921 | 0.940 | 0.899 | 0.944 | 0.932 | 1.055 |
| Proportion Female | 0.62 | 0.64 | 0.68 | 0.80 | 0.68 | 0.75 | 0.71 | 0.84 | 0.73 | 0.77 | 0.48 |
| Proportion Male | 0.38 | 0.36 | 0.32 | 0.20 | 0.32 | 0.25 | 0.29 | 0.16 | 0.27 | 0.23 | 0.52 |
| Taxon: sheep | |||||||||||
| Lambda | 1.010 | 0.993 | 0.998 | 0.957 | 1.006 | 0.928 | 0.949 | 0.907 | 0.956 | 0.942 | 1.059 |
| Proportion Female | 0.65 | 0.65 | 0.70 | 0.82 | 0.70 | 0.77 | 0.73 | 0.86 | 0.75 | 0.79 | 0.48 |
| Proportion Male | 0.35 | 0.35 | 0.30 | 0.18 | 0.30 | 0.23 | 0.27 | 0.14 | 0.25 | 0.21 | 0.52 |
Table 5
Bootstrapped (n = 1000) means of predicted herd growth rate, λ (λboot), proportions of female and male sheep and goat in herds when survivorship probabilities derived from empirical culling profiles for Early (EN) and Middle (MN) Neolithic sites are parameterized as offtake rates in the Lefkovitch population projection matrix (Lefkovitch 1965). Values shown here reflect 200 years of variation in mortality and fertility rates given in Table 1.
| BENKOVAC-BARICE MN | GRADUŠA LOKVE—ISLAM GRČKI MN | SMILČIĆ EN | SMILČIĆ MN | ZEMUNIK DONJI MN | |
|---|---|---|---|---|---|
| Taxon: goat | |||||
| λboot | 1.007 | 0.979 | 0.964 | 1.000 | 0.979 |
| Proportion Female | 0.59 | 0.68 | 0.65 | 0.59 | 0.62 |
| Proportion Male | 0.41 | 0.32 | 0.35 | 0.41 | 0.38 |
| Taxon: sheep | |||||
| λboot | 1.016 | 0.991 | 0.975 | 1.011 | 0.989 |
| Proportion Female | 0.61 | 0.70 | 0.67 | 0.61 | 0.64 |
| Proportion Male | 0.39 | 0.30 | 0.33 | 0.39 | 0.36 |
Table 6
Results of Levene’s test for equality of variances comparing predicted annual herd growth rate (λ) under each strategy’s unadjusted culling rates with actual annual herd multiplication rates (m) after culling rates were optimized. When λ < λmin (0.936) offtake rates were decreased; when λ > λmax (1.005) offtake rates were increased to achieve m = 1, such that herd size remained unchanged from one year to the next. Significant (p < 0.05) results in bold.
| STRATEGY | GOATS | SHEEP | ||||
|---|---|---|---|---|---|---|
| F | p-VALUE | λboot | F | p-VALUE | λboot | |
| Redding (1981) | ||||||
| Energy | 3.983 | 0.047 | 0.992 | 6.130 | 0.014 | 1.010 |
| Security | 4.415 | 0.036 | 0.985 | 5.456 | 0.020 | 0.993 |
| Payne (1973) | ||||||
| Meat | 4.273 | 0.039 | 0.989 | 6.593 | 0.011 | 0.998 |
| Milk | 5.300 | 0.022 | 0.950 | 2.292 | 0.131 | 0.957 |
| Wool | 5.222 | 0.023 | 0.996 | 4.003 | 0.046 | 1.006 |
| Vigne and Helmer (2007) | ||||||
| Meat A | 2.823 | 0.094 | 0.921 | 0.711 | 0.399 | 0.928 |
| Meat B | 1.355 | 0.245 | 0.940 | 3.241 | 0.073 | 0.949 |
| Milk A | 6.899 | 0.009 | 0.899 | 8.615 | 0.004 | 0.907 |
| Milk B | 0.905 | 0.342 | 0.944 | 2.715 | 0.100 | 0.956 |
| Fleece | 7.363 | 0.007 | 0.932 | 7.368 | 0.007 | 0.942 |
| Neolithic | ||||||
| Benkovac-Barice MN | 3.775 | 0.053 | 1.007 | 7.157 | 0.008 | 1.016 |
| Graduša Lokve—Islam Grčki MN | 2.695 | 0.101 | 0.979 | 5.150 | 0.024 | 0.991 |
| Smilčić EN | 1.665 | 0.198 | 0.964 | 2.520 | 0.113 | 0.975 |
| Smilčić MN | 3.599 | 0.059 | 1.000 | 5.810 | 0.016 | 1.011 |
| Zemunik Donji MN | 2.464 | 0.117 | 0.979 | 2.700 | 0.101 | 0.989 |
| No Offtake | ||||||
| Baseline | 0.000 | 1.000 | 1.055 | 0.000 | 1.000 | 1.059 |

Figure 3
Comparison of bootstrapped herd growth rate (λboot) for sheep and goats. Vertical ranges reflect 95% confidence intervals. Horizontal line is shown at λ = 1, where population size is stable.

Figure 4
Initial age structure (combined males and females) of herds of 150 goats and 150 sheep for each culling strategy examined. EN = Early Neolithic, MN = Middle Neolithic.

Figure 5
Goat and sheep population projections over a 200-year period for each culling strategy. The Baseline model projects unconstrained herd population dynamics. Herd sizes start at 150 sheep and 150 goats.

Figure 6
Results of 100 replicated simulations of herd size changes over 200 years for each culling strategy. Initial herd size is set to 150 each for sheep and goats. Each simulation produced a unique set of inter-annual fertility and mortality rates separately for goats and sheep. Replications are consistent with the predictions associated with λboot: growth when λboot > 1, decline when λboot < 1, and steady state is reached when λboot = 1.

Figure 7
Reprojected population dynamics with optimized offtake rates. Optimization prolongs herd survival for all strategies where λboot < 1 and causes a population decline for strategies where λboot > 1.

Figure 8
Mean herd size over the 200-year simulation of goat and sheep population size changes with optimized offtake rates. Mean number of sheep herd sizes are larger than goats for all strategies except for Payne’s (1973) Milk and Vigne and Helmer’s (2007) Meat A strategies.

Figure 9
10-year moving average of annual goat population growth rate with unadjusted offtake rates (λ) versus actual population growth rate after optimization of offtake rates (m). Triangles and inverted triangles show high and low points, respectively. Solid and dashed lines refer to lower (λmin = 0.936) and upper (λmax = 1.005) thresholds used to initiate optimization of culling rates. The degree to which m and λ overlap is an indication of the frequency by which culling rates were adjusted. For example, no overlap is shown for Milk A since λ was almost always less than λmin.

Figure 10
10-year moving average of annual sheep population growth rate with unadjusted offtake rates (λ) versus actual population growth rate after optimization of offtake rates (m). Symbology as in Figure 9.
Table 7
Results of sensitivity analysis reporting Sobol’s first order effects of female offtake on average herd size and interannual variation (sd).
| TAXON | MEAT | MILK | WOOL | |||
|---|---|---|---|---|---|---|
| SIZE | sd | SIZE | sd | SIZE | sd | |
| goat | 0.16 | 0.87 | 0.12 | 0.72 | 0.15 | 0.87 |
| sheep | 0.21 | 0.94 | 0.22 | 0.97 | 0.18 | 0.88 |

Figure 11
First-order indices of sensitivity of mean herd size (size) and interannual variation (sd) to female offtake rates, λmin, and λmax.

Figure 12
Total-effect indices of sensitivity of mean herd size and interannual variation (sd) to female offtake rates, λmin, and λmax.
