In this research, was used the IBM type Agent Based Modeling (ABM). The algorithms used (Supplementary Material I, II and III) were developed within the platform Matlab® version R2015a50. The model description is structured according to the protocol ODD (Overview, Design concepts, Details) updated suggested by Grimm et al. (2010).
Purpose
The purpose of the models is to assess the effect of varying intrinsic parameters on the population viability of an endemic fish species. Although endemism and species unity are not a limiting feature for replicating models. Thus, the developed algorithms indicate how (i) the ratio between the birth and mortality rates (Supplementary II – MetaZebra01) and (ii) the level of specialization and competition (Supplementary III – MetaZebra02) of the specimens affect population viability.
Entities and state variables
The model has only one entity, specimens. The specimens are representatives of a single species of fish distributed in three populations of different sizes. Each specimen was characterized by a set of specific parameters based on information from literature and breeders of the species. Basically, we used characteristics that influence population dynamics, including the following state variables: longevity, age of sexual maturity, annual reproduction number, instant birth rate (b0), instant mortality rate (d0), interference of each individual in population growth (b1), interference of each individual in population mortality (d1), time (t), radius, richness (S), tolerance (tol) and population size (N). We consider the following assumptions for determining the values of our variables:
Longevity. Species of the same family (e.g. Ancistrus ssp.) can reach more than 15 years of age 51–53. According to aquarists, longevity in captivity is at least 15 years for the species under study. However, probably in a natural environment, specimens have a shorter life span when compared to individuals in good health in captivity54. Therefore, in the model, we estimate that specimens can reach 12 years of life in nature.
Age of sexual maturity. According to breeders who produce the species in captivity, individuals reach sexual maturity at the age of three.
Annual reproduction rate. To determine b0, we need to know the number of individuals who enter the cohort each year by birth. Because the male copulates with more than one female55, we only consider females in the calculation. We used 50% of the specimens of each hypothetical population, since the parameter is more related to the fertility of females, as the number of individuals depends more on that sex, so the sex ratio considered was 1:1. Each female spawns an average of 14 eggs per spawn48,55. In captivity, multiple spawning are observed throughout the year. In a natural environment, reproduction can occur at any time of the year, but two reproductive peaks were observed annually56,57. We estimate that 95% of the females in each cohort are able to reproduce. Most females spawn twice a year, possibly larger females with better nutritional performance, are more apt for a greater number of annual reproduction58. We consider the following rates: 35% of females reproduce only once / year, 45% reproduce twice/ year and only 15% of females reproduce three times / year.
Mortality rate. Regarding mortality, we arbitrarily define that in a natural environment despite parental care, the mortality of individuals under one year of age (juvenile and egg) is high, around 50%, due to predation and competition for hide. Additionally, adult females are assigned a 25% mortality rate, while adult males have 20%. We consider that the behavior of the male to stay hidden and protecting it’s crevice (Ramon, 2011; Gonçalves, 2011) provides less vulnerability to males and consequently a lower mortality rate than females in the population.
Intrinsic birth rate (b0). According to the assumptions of the hypothetical natural population size and annual reproduction rate, we obtained the result of b0= 0,6. We used only the females of each population (50% of the individuals) and the rate value was obtained through a difference equation presented in more detail in section 2.6 (sub-models).
Intrinsic mortality rate (d0). Following the assumptions of the hypothetical natural population size and annual mortality rate it was estimated as d0= 0,3. We used only the females of each population (50% of the individuals) and the rate value was obtained through a difference equation presented in more detail in section 2.6 (sub-models).
Interference of each individual in the growth (b1) and mortality (d1) of the population. We did not consider the effect of b1 (interference of each individual on population growth) and d1 (interference of each individual on population mortality) in the models.
Time (t). Time was measured in years in models with t = 0 as a starting point, ending in 1,000 years.
Radius. The radius rate varies from 0 to 1, being 0 when there is no competition and 1 when everyone is competing with each other. There are no studies that say how competitive the species is. So, we determine an average value (0.5).
Richness (S). As the study deals with the population of just one species, we considered S = 1.
Tolerance (tol). There is also no information on the species tolerance, therefore, we determined an average value (0.5).
Initial population size (N). In both modules of model execution, we considered n the initial size of the population with 100,000 individuals.
Process overview and scheduling
The processes of the models promote the simulation of the dynamics of individuals within the population in an environment without anthropic effects and interspecific interaction (Figure 1; Supplementary I – PopZebra). The process begins with the entry of a cohort with an initial number of individuals (n = 10.000) in the population. Gradually, individuals are assigned to the Optimal level range (OLR) randomly determined, or to the Maximum or Minimum tolerance level range (LRMaxMin) ranging from zero to one, according to the model. Specimens included in the OLR are aged, go through the update of age and later young individuals go through the process of sexual maturity until they reach the age of sexual maturity (adult individuals), the individuals enter the reproduction process, causing the origin of a new cohort by birth.
As for the individuals that enter FNMaxMin, they are destined to the competition processes for resources (territorialism and / or food). Randomly, some individuals are classified as survivors and enter the aging process and the consecutive ones mentioned above, while the rest are removed from the model by the mortality process. Individuals who reach the age of longevity are also removed by the model through the process of natural mortality (by age).
This dynamic is generated through the PopZebra program (supplementary material I). However, to meet our objectives, two metaprograms were developed. MetaZebra 01 (supplementary material II) for objective one and MetaZebra 02 (supplementary material III) for objective two. The values of the variables used in the algorithms (shown in Table 3) are based in knowledge gathered in specialized hobbyists’ magazines, hobbyists and fishermen personal communications, as well as experiments carried out in captive breeding program in the laboratory. The MetaZebra 01 algorithm was built to create combinations of birth and mortality, thus varying b0 and d0 from zero to one in 0.05 intervals. Considering H. zebra as r or k strategist (depending on the mortality rate). This procedure will build an interface of values where the x-axis will be the entire variation of the mortality rate while the y-axis will be the entire variation of the birth rate and the 441 cells (21 birth values multiplied by 21 mortality values) will represent all possible combinations between these two rates. Thus, indicating the effect of the ratio between birth and mortality rates on population viability.
Table 3
Condition variables used in the algorithms to generate the models. Randomized values were determined by “Minimum Value: Interval: Maximum value”.
Variables | PopZebra | MetaZebra 01 | MetaZebra 02 |
Maturidade sexual (year) | 3 | - | - |
Longevidade (year) | 12 | - | - |
b0 | - | 0:0.05:1 | 0.6 |
b1 | 0 | - | - |
d0 | - | - | 0.3 |
d1 | 0 | 0:0.05:1 | - |
t (year) | 1,000 | - | - |
raio | - | 0.5 | 0:0.05:1 |
tol | - | 0.5 | 0:0.05:1 |
sp ou S | 1 | - | - |
n | 100,000 | - | - |
As for the MetaZebra 02 algorithm, it was created to vary the tolerance and the radius of competition, thereafter the tol and the radius vary from zero to one in 0.05 intervals. Considering the species as generalist or specialist (depending on the tolerance rate) and little or very competitive (depending on the competition radius rate) in the face of changes. This procedure will generate an interface of values, where the x axis will be the entire variation of the competition radius rate and the y axis will be the entire variation of the tolerance rate and, the 441 cells (21 tolerance values multiplied by 21 radius values competition) will represent all possible combinations between these two rates. This way, it is possible to compare the effect of the species' level of specialization (tolerance) and relate it to the influence of intraspecific competition on the persistence of H. zebra. Each algorithm generated 441 combinations of values (scenarios) that represent the population's probability of survival over an interval of 1,000 years. This probability was calculated from five replicates of each combination.
Design concepts
Basic principles. Population dynamics are maintained with the constant influence of several factors. Growth is determined by the number of entities (individuals) that enter (birth and migration) and leave (mortality and emigration) the population59. The population grows exponentially until it is controlled by the amount of resources available in the environment, called support capacity60. In addition, ecological factors within the habitat will influence the interaction between individuals (intraspecific interaction) and with the environment (Law of Tolerance), which depending on the level of specialization of the species will determine the growth and survival of this set of organisms 61–63. We assumed, in the model under study, the population as being a closed one (entry of individuals by birth and exit by mortality), since it is an endemic species with no migratory characteristics, restricted to about 170 km of the river section.
Emergency. Population dynamics included an emerging result of behavior and interactions between individuals and their habitat.
Interaction. Only internal interactions are recorded in the process. Intraspecific competition is about organisms of the same species that start competing for resources when demand is greater than what the environment can offer, causing a reduction in population growth due to mortality 64,65. Another interaction that will influence the survival of individuals is the level of specificity in terms of habitat (general and specialist species). Species with a more generalized niche tend to occupy large geographical areas and are vulnerable to factors that reach large scales, while specialist species have a restricted niche, occupying small geographical areas, their populations are small and more vulnerable to extinction 63,66,67, usually endemic and rare species are specialists.
Iniciation
Initially, the condition variables b0, b1, d0, d1, radius, time, n, S and tol generate the behavior curve of the population of the single species under study, executed through the IND function, creating the intrinsic values of the species. Then, the characteristics of the individuals are inserted, using information from the following variables: longevity, age of sexual maturity, tolerance value (ranging from tol = 1 to 0 with 0.05 intervals) and random optimal values (random numbers). Optimal values will not influence the model, as they are irrelevant when S = 1. This way, we created initial conditions to generate the first cohort with n = 10,000 individuals aged = zero (year of age). Each individual is created based on the variation of the optimum and tolerance.
Submodels
Intrinsic birth ( b0 ) and mortality ( d0) rates. For one of our models, we had to define b0 and d0 of the estimated natural populations. b0 (Equation 1) and d0 (Equation 2) were calculated using the following expressions59:
\(b=\frac{B}{N.t}\) Equation 1
where B represents the number of births, and N.t represents the size of the current population.
\(d=\frac{D}{N.t}\) Equation 2
Where, D represents the number of deaths, and N.t represents the size of the current population. We considered only the females (50% of the value) of each population in the calculations.
Both b0 (0.6) and d0 (0.3) showed the same result in the three populations.
Aging of individuals and competition. Each cohort is aged by adding one year. Effective number (Ne) is calculated by the product of the radius and the population size (sum of all living individuals in the base year).
Birth. From the third year of life, when the age of sexual maturity begins, individuals start reproduction and births are calculated by individuals, following equation 359:
B=b0+b1*Ne Equation 3
Where, b0 represents the intrinsic birth rate, b1 when each individual interferes with population growth and Ne the effective number.
Mortality due to the number of individuals. Depending on the radius and population growth, individuals are removed from the population, following the following equation 4:
D=d0+d1*Ne Equation 4
Where, d0 represents the intrinsic mortality rate, d1 when each individual interferes in the mortality of the population and Ne the effective number.
Mortality by age. When the surviving individuals of each cohort reach 12 years of age, they are eliminated from the population by the aging mortality process.
Delimitation of natural populations and viability study
As for the third objective, to infer about the viability status of natural populations, we defined the existence of three populations of H. zebra considering the proximity of the known occurrence points and possible natural barriers (Figure 2), since there are no studies that define the size of existing natural populations. Taking into account that we do not know the location of the geographical barriers that separate the populations, we defined only three natural populations: Gorgulho da Rita (from the city of Altamira to Volta Grande), Volta Grande (from the downstream of Altamira to the upstream of Vila Belo Monte) and Belo Monte (from Vila Belo Monte to Vitória do Xingu; Figure 2).
Population sizes were defined based on literature and informal reports by fishermen and personal observations. This informal estimate was supported by the rate of offered smuggled specimens on the internet. According to some information, in 2017, around 10,000 specimens of H. zebra were removed, from nature, monthly in the dry season (July to November). It was assumed that in the high-water period, 40% (4,000 individuals / month) of this quantity was illegally collected. We assume that the number of individuals taken from nature represents 10% of the population, generating a total size of 7,800,000 distributed among three populations. The abundance of each population was estimated based on the abundance of the capture samples of Gonçalves (2011). Given the total size of the species' metapopulation, we consider that 15% (1,170,000 individuals) belong to Gorgulho da Rita, 60% (4,680,000) are in Volta Grande and 25% (1,950,000) in Belo Monte. So, the study of the viability of the three populations was carried out through the interface that combines birth with mortality rates, using the MetaZebra 01 algorithm. We used the values that can be found in equations 1 and 2 (see Submodels) of intrinsic birth rate (b0 = 0.6) and intrinsic mortality rate (d0 = 0.3) for the three natural populations and related their viability status.