Bayesian Parameter Estimation of Cannibalism Dynamics in Larval Atractosteus tropicus Under Different Feeding Conditions

1. Introduction

Cannibalism in fish is an intraspecific interaction that affects survival, size structure, and population. Trophic dynamics, and their importance varies among species and ecological contexts [1]. In cultured fish populations, where space is limited, intracohort cannibalism is often intensified by rapidly growing individuals and by the morphological constraints associated with prey ingestion [2]. For example, in sea bass hatcheries, increased size heterogeneity and feeding frequency have been associated with a higher incidence of cannibalistic events and, consequently, with substantial reductions in survival [3]. In larviculture, cannibalism is now recognized as being more prevalent than previously assumed and may be regarded as an alternative trophic strategy under resource-limited conditions that contributes to species persistence [4]. Under natural conditions, several studies suggest that cannibalism may be rare in certain groups, such as livebearing fishes, but may increase when competition for resources intensifies, particularly at high population densities [5].

Cannibalism may also arise from reproductive conflicts in species exhibiting parental care, in which the consumption of offspring may result from parental competition rather than from an energy deficit [6]. Recent studies have emphasized that population-level and environmental factors may influence the occurrence of cannibalism in a nonlinear way [7], and that increasing the structural complexity of the habitat may promote cannibalism, the coexistence of individuals of different ages, and population growth [8]. In habitats invaded by superpredators, cannibalism may serve as a mechanism that mitigates resource scarcity and contributes to the self-sustainability of invasive populations in human-modified ecosystems [9].

Cannibalism has been investigated mathematically through several models developed to describe population dynamics that are not captured by classical predator–prey models. An approach is to incorporate age structure by dividing the predator population into juvenile and adult classes and introducing cannibalism into the original predator–prey model [10]. Some models have also been proposed in which cannibalism occurs simultaneously in both prey and predator populations within a ratio-dependent Holling–Tanner framework. These models show that, although cannibalism restricted to only one species may induce instability, its simultaneous occurrence in both species may promote their coexistence [11]. Another approach involves Leslie–Gower-type models that incorporate a prey cannibalism term. Such models show that cannibalism does not alter the local stability of the coexistence equilibrium and that conditions for global stability can be established [12].

More complex ecological scenarios have also been simulated, including cannibalism in trophic webs composed of phytoplankton, zooplankton, and fish under proportional harvesting. In this setting, local and global stability are analyzed, and the existence of a Hopf bifurcation with respect to a cannibalism-related parameter is demonstrated [13]. When cannibalism is restricted to specific life stages, such as the juvenile stage, conditions can be obtained under which cannibalism stabilizes the population, as well as conditions under which stabilization cannot be achieved [14]. These findings reinforce the idea that the effect of cannibalism critically depends on how population structure, energy gain, and interaction mechanisms are modeled.

In this work, we perform a Bayesian parameter estimation for a mathematical model based on a system of ordinary differential equations describing the dynamics of a larval population of Atractosteus tropicus reared in a pond. The model divides the population into three subclasses: vulnerable individuals (x), cannibalistic individuals (y), and non-vulnerable individuals (z) [15]. The vulnerable population consists of smaller larvae. The cannibalistic population comprises individuals that are larger than the vulnerable larvae or possess morphological characteristics that enable them to consume conspecifics. The non-vulnerable population consists of larvae that have attained sufficient morphological development to avoid being consumed and that do not possess the characteristics required to consume individuals of their own species. The model is based on the following assumptions:

1.

Initially, the entire larval population belongs to the vulnerable class.

2.

The interaction between the vulnerable and cannibalistic populations is modeled as a predator–prey relationship using a Beddington–DeAngelis functional response.

3.

The quality and acceptance of the supplied food $\left({\alpha }_1\right)$ influence the occurrence of cannibalism.

4.

A rate $\left(m\right)$ for monitoring and isolation of larvae exhibiting cannibalistic behavior is considered.

5.

A proportion of the larval population becomes cannibalistic at rate $\left(\beta \right)$ for a certain period of the time.

6.

All population classes are subject to a mortality rate $\left(\mu \right)$.

The model is specifically given by

$$\begin{array}{cc}\hfill \dot{x}& =-{\displaystyle \frac{(1-{\alpha }_1)\beta xy}{1+x+ay}}-(1-{\alpha }_1)(1-m)\beta x-{\rho }_{1}x-\mu x,\hfill \\ \hfill \dot{y}& =(1-{\alpha }_1)(1-m)\beta x-{\alpha }_1\gamma y^2-{\rho }_{2}y-\mu y,\hfill \\ \hfill \dot{z}& ={\alpha }_{1}m(n-x-z)({\rho }_{1}x+{\rho }_{2}y)-\mu z,\hfill \end{array}$$

(1)

where

• m is the parameter associated with handling time and represents the proportion of each day devoted to monitoring and isolating cannibalistic individuals, with $m\in (0,1)$.
• a quantifies interference among cannibalistic individuals.
• n represents the initial vulnerable population.
• ${\alpha }_1$ quantifies the quality and acceptance of the supplied food, with ${\alpha }_1\in (0,1)$.
• $\mu$ represents the mortality rate of Atractosteus tropicus larvae.
• $\beta$ is the rate at which vulnerable individuals acquire cannibalistic characteristics.
• ${\rho }_1$ and ${\rho }_2$ are the rates at which individuals in the vulnerable and cannibalistic populations, respectively, transition to the non-vulnerable class.
• $\gamma$ is the cannibalism rate within the cannibalistic population.
• The function ${\displaystyle \frac{(1-{\alpha }_1)\beta x}{1+x+ay}}$ describes the attack rate of cannibalistic individuals, regarded as predators, on the vulnerable population, regarded as prey. This expression incorporates the effects of predator satiation and interference.
• The factor $n-x-z$, which appears in the growth term of population z, determines the maximum growth threshold of the non-vulnerable population.

Furthermore, we assume that all parameters of system (1) are positive, and we are interested in the behavior of its solutions during the first seven days following the onset of cannibalism.

2. Parameter Estimation

To carry out the numerical estimation of the mathematical model parameters, the Montecarlo method is used. To achieve this, the mathematical model is solved numerically over multiple iterations, allowing for the subsequent Bayesian estimation of the parameters.

In this work the mathematical system is solved by using the integrator implemented in the LSODA (Livemore Solvers for Ordinary Differential Equation with Automatic Method Switching) library [16]. The library LSODA was written in FORTRAN in the Lawrence Livermore National Laboratory. This algorithm efficiently handles both non-stiff and stiff problems by automatically switching between the Adams-Moulton method and backward differentiation formulas (BDF). This ensures numerical stability when evaluating regions where the population dynamics change rapidly.

2.1. Likelihood Function

Let us suppose the following data set $(t_1,y_{o1}),(t_2,y_{o2}),\dots (t_n,y_{on})$ and the solutions given by the mathematical model: $(t_1,y_{m1}),(t_2,y_{m2}),\dots (t_n,y_{mn})$, where $y_{mi}$ is the $ith$ component of the solution computed for the system

$${\displaystyle \frac{d\mathbf{y}}{dt}}=f(t,\mathbf{y},\theta ).$$

Where $\theta$ is the set of parameters in the mathematical model. In practice there exist differences between the observed data and the computed solution, which reads

$$y_{oi}=y_m(t_i,\theta )+{\epsilon }_i,$$

(2)

where ${\epsilon }_i$ is an error. These errors come from the simplifications the in model, data measurements, and biological effects, among them.

Let us suppose

$${\epsilon }_i\sim N(0,{\sigma }^2),$$

that means, the probability of observing $y_{oi}$ given the model following probability density

$$P\left(y_{oi}|\theta \right)={\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}{e}^{-{\displaystyle \frac{{(y_{oi}-y_m(t_i,\theta ))}^2}{2{\sigma }^2}}}.$$

(3)

The likelihood is the joint probability of observing every all data points given a parameter set $\theta$. It is supposed that errors are independent and they follow the same probability distributions, then

$$\mathcal{L}\left(\theta \right|D)={\displaystyle \prod _{i=1}^n}P\left(y_i\right|\theta ),$$

$$\mathcal{L}\left(\theta \right|D)={\displaystyle \prod _{i=1}^n}\left[{\displaystyle \frac{1}{\sqrt{2\pi }\sigma }}{e}^{-{\displaystyle \frac{{(y_{oi}-y_m(t_i,\theta ))}^2}{2{\sigma }^2}}}\right].$$

(4)

Therefore, the log likelihood read as

$$ln\left(\mathcal{L}\left(\theta \right|D)\right)=-nln\left(\sigma \sqrt{2\pi }\right)+{\displaystyle \sum _{i=1}^n}-{\displaystyle \frac{{(y_{oi}-y_{mi})}^2}{2{\sigma }^2}}.$$

(5)

According to [20], if $\sigma$ is fixed, then the term $-nln\left(\sigma \sqrt{2\pi }\right)$ is not taken into account for the compute the maximum value of the log likelihood function (5). Hence, the log-likelihood becomes in

$$ln\left(\mathcal{L}\left(\theta \right|D)\right)={\displaystyle \sum _{i=1}^n}-{\displaystyle \frac{{(y_{oi}-y_{mi})}^2}{2{\sigma }^2}}.$$

(6)

2.2. Bayes Theorem

The Bayes Theorem read as

$${\displaystyle \mathbf{P}\left(\theta \right|D)=\frac{\mathbf{P}\left(D\right|\theta \left)\mathbf{P}\right(\theta )}{\mathbf{P}\left(D\right)},}$$

(7)

where

• ${\displaystyle \mathbf{P}\left(\theta \right|D)}$ is the posterior probability function.
• ${\displaystyle \mathbf{P}\left(D\right|\theta )}$ is the likelihood function.
• ${\displaystyle \mathbf{P}\left(\theta \right)}$ is the prior probability function.
• ${\displaystyle \mathbf{P}\left(D\right)}$ is the data probability.

Applying the Bayes Theorem and taking into account (6), we can get

$${\displaystyle ln\left(\mathbf{P}\right(\theta \left|D\right))=ln(\mathcal{L}\left(\theta \right|D\left)\right)+ln\left(\mathbf{P}\right(\theta \left)\right)-ln\left(\mathbf{P}\right(D\left)\right).}$$

(8)

By supposing the independence for the parameter in the mathematical model the joint probability function for $\theta$, reads

$${\displaystyle \mathbf{P}\left(\theta \right)={\displaystyle \prod _{i=1}^N}{f}_j\left({\theta }_j\right|{\eta }_j),}$$

(9)

being

1.

N is the total of parameters for the model.

2.

${\theta }_j$ is the j-th parameter for the model.

3.

$f_j$ in the density probability function for ${\theta }_j$.

4.

${\eta }_j$ is the set of parameters for the density probability function $(\mu ,\sigma ,[a,b\left]\right)$.

Therefore

$$ln\left(\mathbf{P}\left(\theta \right)\right)={\displaystyle \sum _{j=1}^N}ln\left(f_j\left({\theta }_j\right|{\eta }_j)\right).$$

(10)

Following the works [21,22], the acceptance probability, in the Metropolis-Hastings algorithm, depends on the ratio of the posterior densities. Since $\mathbf{P}\left(D\right)$ is independent of the data, it is canceled in the Metropolis ratio. We can get

$$ln\left(\mathbf{P}\left(\theta \right|D)\right)\propto -{\displaystyle \frac{1}{2}}{\displaystyle \sum _{i=1}^n}{\left({\displaystyle \frac{y_{oi}-y_{mi}}{\sigma }}\right)}^2+{\displaystyle \sum _{j=1}^N}ln\left(f_j\left({\theta }_j\right|{\eta }_j)\right).$$

(11)

2.3. Metropolis - Hastings Method

This method look for to generate a sequence of samples $\{{\theta }^{\left(0\right)},{\theta }^{\left(1\right)},\dots {\theta }^{\left(N\right)}\}$ such that it converges to the stationary distribution $\pi \left(\theta \right|D)$. This stationary distribution will be the posterior distribution.

Let ${\theta }^{\left(0\right)},{\mathbb{R}}^k$ such that

$$\mathbf{P}\left({\theta }^{\left(0\right)}\right|D)>0,$$

where k is the total number of parameters to be determined. For each iteration t, a new state ${\theta }^{(*)}$ is generated from the state ${\theta }^{\left(t\right)}$ by using a proposal distribution

$$q\left({\theta }^{(*)}\right|{\theta }^{\left(t\right)})=N({\theta }^{(*)};{\theta }^{\left(t\right)},{\Sigma }_{\delta }).$$

(12)

Where ${\Sigma }_{\delta }$ is a square matrix of size k, such that

$${\Sigma }_{{\delta }_{ij}}=\left\{\begin{array}{ccc}\hfill {\delta }_i^2& \mathrm{if}& i=j,\hfill \\ \hfill 0& & \mathrm{else}.\hfill \end{array}\right.Fori<k\mathrm{and}j<k.$$

Following [20], by using (12), the new state ${\theta }^{(*)}$ is read as

$${\theta }_j^{(*)}={\theta }_j^{\left(t\right)}+\epsilon ,$$

(13)

for $\epsilon \sim N(0,{\delta }_i^2)$ being $\delta$ the step size. For a target distribution $\pi \left(\theta \right)$, the optimal acceptance rate is around of 23.43% [23,24].

2.4. Acceptance Probability

The evaluation of the posterior distribution $\pi \left({\theta }^{(*)}\right)$ for every proposed candidate requires the value of the likelihood function (6). To compute it, the LSODA integrator is called to solve the ordinary differential equations system. Then, the computed solution $y_m(t_i,{\theta }^{(*)})$ is employed to evaluate (6), that solution is also to evaluate the log-prior function (8).

The acceptance probability $\alpha$ is defined by using the Hasting ratio as

$$\alpha \left({\theta }^{\left(t\right)},{\theta }^{(*)}\right)=min\left(1,{\displaystyle \frac{\pi \left({\theta }^{(*)}\right)q\left({\theta }^{\left(t\right)}\right|{\theta }^{(*)})}{\pi \left({\theta }^{\left(t\right)}\right)q\left({\theta }^{(*)}\right)|{\theta }^{\left(t\right)}}}\right)$$

(14)

Since the normal distribution is symmetric, ${\displaystyle q\left({\theta }^{\left(t\right)}\right|{\theta }^{(*)})=q\left({\theta }^{(*)}\right|{\theta }^{\left(t\right)})}$.Then, we get

$$\alpha \left({\theta }^{\left(t\right)},{\theta }^{(*)}\right)=min\left(1,{\displaystyle \frac{\pi \left({\theta }^{(*)}\right)}{\pi \left({\theta }^{\left(t\right)}\right)}}\right).$$

(15)

We can define the log ratio as

$$r_{ln}=ln\left(\pi \left({\theta }^{(*)}\right)\right)-ln\left(\pi \left({\theta }^{\left(t\right)}\right)\right),$$

(16)

In order to guarantee the balance detailed condition, let u be a random number such as $u\sim U(0,1)$

• If $\pi \left({\theta }^{(*)}\right)>\pi \left({\theta }^{\left(t\right)}\right)$, then $r_{ln}>0$. Therefore the acceptance probability is equal to one and ${\theta }^{(*)}$ is accepted.
• If $\pi \left({\theta }^{(*)}\right)<\pi \left({\theta }^{\left(t\right)}\right)$, then $r_{ln}<0$. Hence, ${\theta }^{(*)}$ does not improve upon ${\theta }^{\left(t\right)}$.
  • If $ln\left(u\right)<r_{ln}$ then ${\theta }^{(*)}$ is accepted.
  • If $ln\left(u\right)>r_{ln}$ then ${\theta }^{(*)}$ is not accepted.
  This condition allows the scheme to escape local maximums, and it explores the whole distribution.

Then

$${\theta }^{(t+1)}=\left\{\begin{array}{ccc}{\theta }^{(*)}& if& {\theta }^{(*)}\mathrm{is}\mathrm{accepted},\\ {\theta }^{\left(t\right)}& if& {\theta }^{(*)}\mathrm{is}\mathrm{rejected}.\end{array}\right.$$

(17)

3. Numerical Results

In this section we develop the parameter estimation for two sets of data, these sets of data were collected from [25]. In that work, they have studied the interaction of larvae in two cases: when the larvae were fed and when they were not. Those two cases were studied to measure if cannibalism appears as a result of survival processes or whether it is just an instinctive behavior

3.1. Case 1: Larvae Is Fed

In this first case, we simulated the scenario where the larvae are fed. The obtained experimental data are shown in Table 1.

To carry out the simulation, we developed Python code that implements the Metropolis-Hastings method for estimating the parameters of the mathematical model.

The acceptance for this experiment was 25.3%, this means that the method has explored all the possibles values for the parameters. The numerical approximated parameters are shown in Table 2. The information in the named column initial corresponds to the obtained values from the reference experiments, the prior information.

In Figure 1 we can see the obtained solution for every data set in Table 2. Notice that the obtained solution by using every estimated dataset parameter adjusts in good agreement with the experimental data. The differences between these solutions are minimal, but the one obtained using the mean is the best. Furthermore, we can notice how the solution improves from the prior information to the posterior information. In this case, we have proposed a uniform distribution centered on the prior information for every parameter in the mathematical model. Comparisons between the prior distribution and posterior distribution are shown in Figure 2. Also, the posterior distribution for each parameter after Bayesian inference is shown in Table 3, and these distributions are depicted in Figure 3.

To investigate the relationships among the parameters estimated for the fed population, we calculated the correlation matrix shown in Figure 4. Most correlations are weak, suggesting that the parameters represent largely independent biological mechanisms within the model for the fed population. Nevertheless, several important relationships can be identified.

The strongest correlation, with a value of $-0.81$, occurs between ${\alpha }_1$ and m, which represent food quality and handling time, respectively. Biologically, this indicates that improving the quality of the supplied food reduces the need to monitor and isolate cannibalistic individuals. Thus, adequate feeding may reduce the occurrence of cannibalistic behavior in the initial population.

A moderate positive correlation of $0.33$ is observed between $\beta$ and m, indicating that an increase in the rate at which vulnerable individuals transition to the cannibalistic class entails a greater handling time. Likewise, the negative correlations of $-0.31$ between ${\alpha }_1$ and ${\rho }_1$ and $-0.28$ between ${\alpha }_1$ and ${\rho }_2$ indicate that better feeding conditions reduce the need for vulnerable and cannibalistic individuals to transition to the non-vulnerable class. This may occur because high-quality food decreases the stress and competition associated with cannibalism.

Moreover, parameters such as the interference coefficient a, the mortality rate $\mu$, and the cannibalism parameter $\gamma$ exhibit correlations close to zero with the remaining parameters. This suggests that these mechanisms operate relatively independently within the population dynamics of the fed system. Overall, the correlation analysis highlights food quality as one of the most important factors in reducing cannibalistic interactions and the handling time required during the culture of Atractosteus tropicus larvae under fed conditions.

3.2. Case 2: Larvae Is Not Fed

In this section the behavior when the larvae are not fed is simulated, the data for this experiment is shown in Table 4.

In this simulation we reach an acceptance rate of 23.8%, and the approximated numerical parameters are shown in Table 5. From the Figure 5, we can observe an important improvement of the numerical solution given by the prior information in comparison with the one given by the posterior information. In this case, the best approximation in the first time is given by the mode, but the best adjustment for the final ones is the mean; therefore, in order to predict the posterior behavior of the larvae population, both of them can be employed. As we have done before, in this case, a uniform distribution centered on the theoretical value is taken into account as prior information. The comparison between the prior distribution and posterior distribution is shown in Figure 6. In Table 6 we can find the best-adjusted marginal distribution for each of the parameters.

Figure 7 presents the correlation analysis among the parameters for the unfed population. Under these conditions, the larval population decreased from $100\%$ on day 8 to only $27.14\%$ on day 14, demonstrating the strong effect of food deprivation on the survival of Atractosteus tropicus larvae. The correlation matrix indicates that, under nutritional stress, the dynamics of the system are primarily determined by the interplay among mortality, the development of surviving individuals, and the emergence of cannibalistic behavior.

The strongest positive correlation, with a value of $0.83$, occurs between the mortality rate $\mu$ and the transition rate of cannibalistic individuals to the non-vulnerable class ${\rho }_2$. Biologically, this relationship suggests that as mortality caused by food deprivation increases, population persistence becomes more dependent on the ability of cannibalistic individuals to develop into less vulnerable stages. This behavior is consistent with the hypothesis that, in the absence of external resources, only individuals with greater competitive advantages are able to survive and continue their development [26]. This relationship also indicates that the emergence of cannibalistic behavior is one of the most important mechanisms contributing to the persistence of a proportion of the population. In addition, a positive correlation of $0.70$ is observed between the rate at which vulnerable individuals become cannibalistic, $\beta$, and the handling time, m, indicating that a greater prevalence of cannibalistic individuals requires more handling time. The estimated values of $\beta$ reported in Table 5 remain high, whereas the parameter m increases considerably relative to its theoretical value. This suggests that a greater prevalence of cannibalism would require increased monitoring and isolation efforts to reduce population losses under culture conditions.

Figure 8. Better distribution for each one of the parameters.

Figure 8. Better distribution for each one of the parameters.

Conversely, the negative correlations between the cannibalism rate $\gamma$ and the transition rate of vulnerable individuals to the non-vulnerable class ${\rho }_1$ ($-0.26$), as well as between interference among cannibalistic individuals a and ${\rho }_2$ ($-0.25$), indicate that increased aggressive interactions hinder the normal development of individuals. In particular, an increase in cannibalism reduces the likelihood that vulnerable individuals will survive, whereas greater interference among cannibalistic individuals decreases the efficiency with which they exploit the resources obtained through cannibalism. Furthermore, Table 5 shows that the parameter ${\alpha }_1$, which represents food quality, is lower than its theoretical value, reflecting the effective absence of nutritional resources during the experiment. These results suggest that when the larval population is not fed, cannibalism ceases to function solely as a mechanism of intraspecific competition and instead becomes a survival strategy that contributes to the persistence of a small proportion of the population. Therefore, the correlation matrix supports the hypothesis that food availability is one of the most important factors regulating the intensity of cannibalism and maintaining the survival of Atractosteus tropicus larvae in culture systems.

4. Conclusions

In this study, Bayesian parameter estimation using the Metropolis-Hastings algorithm was successfully applied to a mathematical model who describes the cannibalistic dynamics of Atractosteus tropicus larvae. The model was evaluated under two distinct nutritional scenarios, fed and unfed, providing deep insights into the mechanisms driving intracohort cannibalism.

The numerical results highlight that food availability fundamentally alters the role of cannibalism within the population:

Under fed conditions: Cannibalism acts primarily as an opportunistic behavior driven by intraspecific competition. The parameter estimation revealed a strong negative correlation $-0.81$ between food quality, ${\alpha }_1$, and the required handling time m. This demonstrates that supplying high-quality food is the most effective strategy to mitigate cannibalistic interactions, consequently reducing the operational effort required to monitor and isolate aggressive individuals.

Under unfed conditions: The dynamics shift drastically, and cannibalism becomes a critical survival strategy rather than mere competition. This is evidenced by the strong positive correlation $0.83$ between the mortality rate $\mu$ and the transition rate of cannibalistic individuals to the non-vulnerable class ${\rho }_2$. In the absence of external nutritional resources, only those individuals that resort to cannibalism acquire the necessary energy to develop into advanced, non vulnerable stages, allowing a small fraction of the population 27.14% to persist.

Furthermore, both scenarios show a positive correlation between the rate at which vulnerable individuals become cannibalistic $\beta$ and the handling time m. This underscores a direct practical application for larviculture: minimizing the emergence of cannibals through adequate feeding protocols directly optimizes farm management by reducing labor-intensive monitoring.

In summary, the mathematical modeling and Bayesian inference applied in this work quantitatively validate that cannibalism in Atractosteus tropicus is a complex, resource-dependent dynamic. Providing high-quality food substantially reduces both cannibalism and handling time, thereby maximizing the overall survival and morphological development of the cultured species