\documentclass[11pt,a4paper]{article}
\usepackage{authblk} % Package to manage author and affiliations

% Language setting
\usepackage[british]{babel}

% Set page size and margins
%\usepackage[a4paper,top=2cm,bottom=2cm,left=2.5cm,right=2.5cm,marginparwidth=1.75cm]{geometry}

\usepackage[utf8]{inputenc}
\usepackage[utf8]{xcolor}
\usepackage{amsmath}
\usepackage{natbib}
\usepackage{apalike}
\usepackage{epsfig}
\usepackage{booktabs} 
\usepackage{setspace}


\renewcommand\Authfont{\scshape}
\renewcommand\Affilfont{\itshape\small}

\setlength{\affilsep}{20em}


\usepackage{lineno}


\title{Evolutionary Dynamics of Oviposition Strategies in Outbreaking Herbivorous Insect Populations: A Stochastic Spatially Explicit PDE Model}

\author[1,2]{Bruzzone, O.A.}
%\cortext[oab]{Bruzzone, O.A.}
%\ead[1]{okktawio@gmail.com}
\author[3]{Rossini, L.}
\author[4]{Aguirre, M.B.}
\author[5]{Crespo, J.}

% Define affiliations
\affil[1]{Fundacion Patagonia Analytics,    Calle 19,    Bariloche,    8400,    Provincia de Río Negro,     Argentina}
\affil[2]{CONICET, CCT-Patagonia Norte, Av. de los Pioneros 2350, Bariloche, 8400, Provincia de Rio Negro, Argentina}
\affil[3]{Service d'Automatique et d'Analyse des Systèmes, Université Libre de Bruxelles, Av. FD Roosevelt 50, cp 165/55, 1050 Brussels, Belgium}
\affil[4]{FUEDEI, Gral. Simón Bolívar 1559, Hurlingham, 1686, Provincia de Buenos Aires, Argentina}
\affil[5]{IEGEBA, (CONICET/UBA), Ciudad Universitaria, 1428, Ciudad Autónoma de Buenos Aires, Argentina}

\date{}  % Remove date

\begin{document}

\maketitle
\clearpage

\begin{abstract}  
Insects show complex dynamics in evading natural enemies. Outbreaking herbivorous insects often spatially aggregate eggs as protection strategy against predators and parasitoids. While this strategy fosters collective defense, it is not common in non-outbreaking species.  

This study introduces a spatially explicit Partial Differential Equations (PDE) model that includes an inheritable egg aggregation strategy. The main aim is to assess its effectiveness as an Evolutionary Stable Strategy (ESS) and its role in outbreaking dynamics. The model simulates the biological and spatial dynamics of the herbivore population from the egg stage to adult stages, while natural enemies target the eggs.

The relationship between egg aggregation and outbreaks was explored through 256 simulations with fixed low ($\sigma=0.1$), medium ($\sigma=1.0$) and high ($\sigma=10$) aggregation, as well as variable inheritable aggregation ($\sigma=0.1$–$10$), under three types of functional response. We tracked changes in herbivore eggs (E), adults (A), and predator/parasitoid (P) populations, considering allele frequency shifts as well.

The results indicate that type I response supports low aggregation, stabilising herbivore populations, while type II response leads to instability, making variable aggregation advantageous post-outbreak. Type III responses show intermediate effects. The lower $\sigma$ spreads the risk, aiding evasion, while the higher $\sigma$ concentrates the risk, assuming predators/parasitoids have consumption limits on outbreak fringes. This promotes local competition, destabilising populations as risk-spreading individuals are outcompeted.  
\end{abstract}  
\clearpage

\linenumbers

\section*{Keywords}
Adaptive evolution, spatial dispersal, functional responses, ecological modelling, population stability

\section{Introduction}

Egg aggregation, the clustering of eggs in masses or in limited favourable areas, is a widespread strategy among herbivorous insects \citep{hanski1987pine,cappuccino1995spatial}. It is especially common among outbreaking species, which often show gregarious behaviour during their immature stages \citep{cappuccino2000oviposition,battisti2023aggregative,wallner1987factors,hunter1991traits}. To date, only 1.6\% of Lepidoptera are known as outbreaking species \citep{hunter1995ecology} but 5-8\% of them aggregate eggs or form groups \citep{stamp1980egg,hebert1983egg}. These species follow a boom-and-bust population cycle, expanding rapidly and then collapsing \citep{may1976bifurcations}.  

Several explanations for egg aggregation exist. According to the \textit{aggregation of risk} hypothesis \citep{chesson1986aggregation} clustering reduces individual predation risk, as searching and consumption efficiency of predators is diluted. This \textit{dilution effect} \citep{gross1993insect} support the hosts by limiting the efficiency of natural e  nemies, particularly when predators or parasitoids target only the peripheral area of a cluster \citep{allen1990influence}. Additionally, aggregation may reduce detection via visual, olfactory, or defensive mechanisms \citep{murali2019dynamic,eisner2000chemical,nishimura2023protective}, and enhance survival via collective behaviours \citep{stamp1980egg}. Despite beneficial aspects, large aggregations can also attract more predators \citep{ioannou2008searching, pearce2006predators}, suggesting a trade-off between aggregation and dispersal.  

This trade-off should align with an evolutionarily stable strategy (ESS) \citep{smith1972,smith1973logic}, where aggregation maximises herbivore fitness by balancing survival benefits and predation risks. While empirical studies offer insights, mathematical models allow an in-depth exploration of these mechanisms. Predator-prey interactions are often modelled through Lotka–Volterra equations \citep{lotka1925elements, volterra1926variazioni} or Nicholson-Bailey models \citep{nicholson1935balance} and focus more on the cyclic population dynamics. Functional response models, instead, describe how predators react to prey density \citep{holling1959some}, usually considering a Type I (linear), Type II (saturating), and Type III (sigmoidal) responses that shape stability \citep{bruzzone2022revisiting}. Type I and II responses destabilise populations \citep{oaten1975functional}, while Type III promotes stability. However, in spatially explicit models, stability can also arise under Type I responses \citep{oaten1975functional}.  

The current literature on modelling oviposition strategies \citep{kershenbaum2012modelling, DOSANJOS2023} lacks the outbreak dynamics, while models describing outbreaks \citep{sharma2015stochasticity, bruzzone2023effect, jepsen2023spatiotemporal, veran2015modeling} discard the oviposition behaviour. The missing link between aggregation and predator responses is a limit from a modelling point of view, as existing models are inadequte to describe this phenomenon properly. Although outbreaking species are more prone to cluster eggs, a direct experimental validation is challenging due to spatial and temporal constraints. Thus, a Lotka-Volterra model incorporating spatially explicit oviposition strategies and predator functional responses is necessary to test this relationship.

An ESS for oviposition is a consequence of natural selection driven by environmental pressures and phenotypic diversity \citep{alger2023evolutionarily, biswas2023evolutionarily}. Evolution supports traits that enhance fitness and shapes genetic variation over time \citep{engen2020phenotypic,allendorf2022}. At the same time, outbreaks reduce genetic diversity through population bottlenecks \citep{kang2023population,grandison2023epidemiology}, potentially affecting the oviposition behaviour. The causal link between genetic diversity loss and egg aggregation is still unexplored, highlighting the need for models that simulate evolutionary responses to predation-induced bottlenecks and outbreaks.  


\subsection{Objectives}  
\noindent  
This study explores, through dynamical modelling, how the functional response of natural enemies affects the evolution of prey egg/host aggregation behaviour and outbreak frequency. The functional response of natural enemies, which defines the consumption rates with respect to local prey/host density, shapes the offspring aggregation trade-offs in space and time. A mathematical analysis of how aggregation affects population dynamics is crucial to improve existing models and to enhance pest control strategies.  

We specifically aim to explore: \textit{i)} if egg aggregation constitutes an evolutionarily stable strategy under different natural enemy functional responses, and \textit{ii)} if it remains advantageous when the efficiency of the enemy declines at high local host densities. To address these questions, we carried out simulations based on a stochastic, spatially explicit Lotka-Volterra model that includes different functional response types. Our investigation has been conducted in two phases. First, we formulated and analysed the stability of spatially explicit models with fixed egg aggregation strategies. Second, we analysed shifts in statistical distribution of aggregation strategies in models that account inheritable variability, testing their evolutionary trajectories under different natural enemy responses.  




\section{Methods}
\noindent
The first step of this study was the formulation of a predator-prey/consumer-resource model based on a stochastic partial differential equation framework that simulates a diffusion process. This model aims at describing the implications of the dilution effect on the efficiency of natural enemies, in the case of large egg aggregates. Our starting point is the eco-evolutionary model of \cite{paez2017eco}, where the evolution of the parameters describes the characteristics of the population, by assumption. Thus, we conducted a simulation in four dimensions: time, genetics, and two spatial dimensions. In our framework, the inheritance of egg aggregation behaviour varies across generations. We also consider the spatially-explicit variation of the Lotka-Volterra predator-prey model proposed by \cite{mclaughlin1991pattern} which, in turn, derives from the Fisher equation in \cite{fisher1937wave}.
A conceptual description of the model is reported in Fig \ref{diagcon}.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.0\textwidth]{SAVE.eps}
  \caption{\label{diagcon} Conceptual Diagram of the Model. The system consists of two populations: herbivores, represented by eggs ($E$) and adults ($H$), and a specific natural enemy that consumes the eggs ($P$). Both adult herbivores ($H$) and natural enemies ($P$) move through the matrix via diffusion, depicted by bell curves in the diagram. Over time, adult herbivores continuously lay eggs, represented by stacked ovals (the taller the stacks, the higher the level of aggregation). Natural enemies reproduce by consuming the eggs.}
\end{figure}


\subsection{Basic Model}
\noindent
A compartmental model \citep{RossiniODE} divides the life cycle of herbivorous species into two stages: eggs ($E$) and adults ($A$). The model also includes a natural enemy that consumes the eggs, denoted by ($P$), and which can generally be a parasitoid or a predator. Mathematically, we obtain the following Lotka-Volterra system of equations:
%
\begin{align}
    \begin{split}
        \frac{dE}{dt} = kA - pEP - eE\\
        \frac{dA}{dt} = eE - d_aA \\
        \frac{dP}{dt} = qpEP - d_pP
    \end{split}
    \label{compartmental}
\end{align}
%
\noindent
where $k$ is the oviposition rate of the herbivores population, $e$ is the egg hatching rate, $p$ is the egg consumption rate of the natural enemy $P$, $q$ is the number of natural enemies that are produced by consuming one herbivore egg, and $d_a$ and $d_p$ are the death rates of the adult herbivore and natural enemy, respectively.


To transform the system (\ref{compartmental}) into a spatially explicit dynamical system we can draw inspiration from the Fisher model \citep{fisher1937wave}. The result is a reaction-diffusion partial differential equation (PDE), where diffusion describes the movement of individuals through space and a reaction term accounts for biological interactions. Reaction-diffusion models have a long tradition in population dynamics \citep{holmes1994}, and have been applied to predator-prey systems since the 1970s (e.g., in \cite{murray1975non}). These models typically generate wave-like patterns in which predators chase prey spatially.

The model ha been extended by including explicit spatial variables, as in \citep{okubo2001diffusion}, in which: \textit{i)} the adults of the herbivores move randomly across a homogeneous area, and \textit{ii)} the eggs remain at the location where they are laid. Mathematically, there is an additional diffusion term into the equations that describe $A$ and $P$:

\begin{subequations}
    \begin{align}
        \frac{\partial E}{\partial t} &= kA - pEP - eE \label{DCE}\\
        \frac{\partial A}{\partial t} &= eE - d_aA + D_a \nabla^2 A \label{DCA}\\
        \frac{\partial P}{\partial t} &= qpEP - d_pP + D_p \nabla^2 P \label{DCP}
    \end{align}
\end{subequations}

\noindent
where $D_a$ and $D_p$ represent the diffusion coefficients for the adult herbivores and the natural enemies, respectively.

\subsection{Model of Aggregation Behaviour in the Herbivores}
\noindent
To better investigate the effect of egg aggregation, let us consider additional assumptions that transform the model from deterministic \eqref{DCE}--\eqref{DCP} into stochastic.



\subsubsection{Conceptual framework for egg aggregation}

\noindent

``Aggregation'', denotes either clustering eggs into a single mass or laying eggs in a reduced area without forming distinct masses. Egg clusters in a single mass aligns with true aggregation. From a landscape perspective, it may have negative spatial auto-correlation similar to blue noise caused by local movement or transportation \citep{de2012blue}, indicating a uniform spatial distribution.
Conversely, even if eggs are evenly distributed across the landscape without forming masses, they may still appear ``aggregated'' if laid preferentially in more favourable areas or near other conspecifics, as in \cite{battisti2023aggregative, goodsman2016aggregation}. To avoid any misunderstanding, this study refers to ``egg aggregate'' as synonymous of ``egg masses'', although they are not necessarily identical phenomena.

We assume a homogeneous distribution of the resources into the space, so that any spatial pattern is due to the movement of herbivores and natural enemies (e.g., the local extinction-colonisation process). This process, as described by \cite{hastings1977spatial, hassell1991spatial, okubo2001diffusion}, leads to the population dynamics characterised by travelling waves \citep{sherratt2001periodic, fu2015wave}. 
Thus, this study refers to ``spatial heterogeneity'' to describe the patterns observed as a result of dynamic interactions.

\subsubsection{Spatial stochastic modelling of egg aggregation behaviour}

\noindent
The model behaves as a spatial stochastic process, where herbivorous insects live and develop in a uniform environment of identical host plants, mathematically represented by a matrix. The egg aggregation behaviour is the stochastic process that indicates the proportion of available eggs laid in a single location within the host plant matrix, while the number of eggs laid in a single cell by a female is the egg mass.

The aggregation behaviour is quantified considering the standard deviation of a gamma probability distribution that describes the egg mass sizes. The standard deviation are collected into a vector that ranges from nearly zero (the minimal variation in egg mass sizes and aggregation) to one (a few large egg masses and many small ones). This vector reflects the variability in the egg aggregation behaviour over the generations, and it is denoted by $\sigma_{\alpha}$.

Egg aggregation is expressed by the number of eggs laid in a single cell with respect to the average number of eggs per cell. Low aggregation occurs when egg masses are nearly identical in size; high aggregation occurs when a few egg masses contain most of the eggs, with many cells having few or no eggs. At each step, the population of herbivores lays a number of eggs described by
%
\begin{equation}
o_s \sim A(x,y,t) \cdot G(\mu_o, \sigma_{\alpha})
\label{govip}
\end{equation}
\noindent
where, $o_s$ is the number of eggs laid, $A(x,y,t)$ is the number of adults in position $(x,y)$ at time $t$, and $\mu_o$ and $\sigma_{\alpha}$ are is the average size and the standard deviation of the egg masses, respectively. The symbol $G$ recalls that the process is gamma-distributed. 

Frequency may shift across the generations of the herbivorous species within a set of genetically-based behavioural patterns. These patterns are the key drivers of cyclic outbreaks and extinctions, and they are subject to natural selection. As a result, less advantageous behaviours decrease while the traits held by the survivors emerge.

It is important to clarify that, according to the objectives of this study, we are not going to consider new genetically-based behaviours. Instead, the simulation starts with a predefined gene pool containing all the possible behaviours. These behaviours are fixed (stereotyped) over the lifetime of the herbivores and are directly connected with specific alleles.

The progeny of the prey or of the hosts, produced by egg aggregation, may inherit the egg aggregation trait, denoted by $\sigma_{\alpha}$, from their parents through a normally-distributed process. Consequently, some offspring may show slight behavioural differences with respect to their parents due to mutations or different gene expression, although remaining within the same genetic pool. As a result, the allele frequency of the aggregation trait, denoted by $\sigma_{\phi}$, will shift in the next generation according to the equation:
%
\begin{equation}
\sigma_{\phi} \sim N(\sigma_{\alpha}, \sigma_h)
\label{hovip}
\end{equation}

\noindent
where $\sigma_h$ is the standard deviation associated with the normal ($N$) egg mass inheritance process (Fig. \ref{hoviposition}).

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.0\textwidth]{Oviposition.eps}
  \caption{\label{hoviposition} Egg Aggregation Inheritance Process. The upper vector represents the aggregation characteristic $\sigma$ of parent adults, while the lower vector represents the same characteristic for their egg offspring in a given matrix cell. The characteristic is inherited with normally distributed variation.}
\end{figure}


\subsubsection{Modelling egg-laying behaviour in herbivores}
\noindent
Let us calculate the probability $p_{\phi\alpha}$ that an adult of the herbivorous species having an aggregation behaviour $\sigma_\alpha$ lays an egg mass in the position $(x,y)$ at time $t$. This probability is given by
%
\begin{equation}
p_{\phi\alpha}(A_{\alpha}(x,y,t)) = 1 - r ^ {A_{\alpha}(x,y,t)}
\label{povip}
\end{equation}

\noindent
where the intrinsic rate of population growth, $r$, is a stochastic process the represents the probability of oviposition when the population size is $A=1$. As the population size increases, the probability of an oviposition event increases and approaches 1.

By merging the equations \ref{govip} and \ref{povip} it is possible to obtain the number of eggs laid by the herbivores over time in each $(x,y)$ cell, namely:
%
\begin{equation}
O(x,y,t) = o_s(x,y,t) (p_{ea}(A_{\alpha}(x,y,t)) > U(0, 1)) 
\label{novip}
\end{equation}
\noindent
where $O(x,y,t)$ is the realised oviposition size at time $t$ and location $(x,y)$ given an estimated oviposition event of size $o_s$ and probability $p_{ea}$; $U(0,1)$ is a uniformly distributed random process from 0 to 1 (inclusive). If $p_{ea} > U(0,1)$, an oviposition with size $o_s$ occurs; otherwise, there is no oviposition and $O$ is zero.

\subsubsection{Modelling predation dynamics and functional responses in egg stocks}
\noindent
We assume that natural enemies consume the eggs following  the functional response in \citet{rogers1972random}, but in a differential equation form as in \citet{rosenbaum2018fitting, bruzzone2022revisiting}. In each position, the stock of egg-prey/host is limited and the eggs in position $x,y$ compose a single stock of prey/host, thus they are assumed to be randomly distributed. The natural enemy will then consume eggs proportionally to their abundance, regardless of the value of $\sigma_{\alpha}$. The functional response equations are thus defined as follows:
%
\begin{subequations}
    \begin{align}
        \text{Holling's type I: } & \frac{dE}{dt} = -aEP & \label{fr1}\\
        \text{Holling's type II: } &\frac{dE}{dt} = -P \frac{aE}{1+ahE}
        \label{fr2}\\
        \text{Holling's type III: } &\frac{dE}{dt} = -P \frac{aE^2}{1+ahE^2}
        \label{fr3}
    \end{align}
    \label{erf}
\end{subequations}

\noindent
Equations (\ref{fr1})--(\ref{fr3}) express the functional response $f_r(E_{x,y,\alpha}P)$ of the natural enemies to the abundance of prey/host in each cell. Equation (\ref{fr1}) is the type I functional response where natural enemies consume prey/host at a constant rate $a$ proportionally to their abundance. Equation (\ref{fr2}) is either the Holling's type-II (or disk) equation in case of unlimited prey/host or a Rogers type II functional response if the stock of prey/host is limited. The additional parameter $h$ is the Holling's handling time. The Equation (\ref{fr3}), instead, is the sigmoidal type III functional response. 

Functional response models are defined through differential equations, so that the prey/host availability is updated at each time step, transforming the overall theoretical framework into the Rogers random predator model \citep{rogers1972random}, as in \cite{rosenbaum2018fitting} and \cite{bruzzone2022revisiting, bruzzone2023new}. In this way, the effect of prey/host depletion is already included in the model.

In the random parasite model in \cite{rogers1972random}, parasitoids cannot distinguish between parasitised and non-parasitised hosts, leading to excessive superparasitism (multiple eggs per host), especially in solitary parasitoids. This results in biologically unrealistic scenarios where some hosts receive an excessive number of parasitoid eggs. Assuming that parasitoids can recognise parasitised hosts, the parasitisation or predation of egg masses is closer to a predator-prey interaction rather than a traditional host-parasitoid relationship. Under these conditions, parasitization directly removes eggs from the prey/host population, akin to predation. As noted by \cite{bruzzone2023new}, using a predator-prey framework is mathematically advantageous and more biologically realistic, particularly when considering host depletion by parasitoids.


\subsection{The complete model}

The equations introduced in the previous sections collapse into a spatially explicit Lotka-Volterra predator-prey model. This model incorporates variable and inheritable oviposition strategies and is more suitable to pursue the objective of this study. By integrating the inherited variable $\sigma_{\alpha}$ into the integro-partial differential form of the first model (Equation \ref{compartmental}), the models are enclosed into a single system:
%
\begin{align}
    \begin{split}
        \frac{\partial E_{x,y,\alpha}}{\partial t} &= k \int_{\alpha = \alpha_n}^{\alpha = \alpha_x} o_s(A_{\alpha}(x,y,t)) d\alpha - f_r(E_{x,y,\alpha}P - hE_{x,y,\alpha}) \\
        \frac{\partial A_{x,y,\alpha}}{\partial t} &= hE_{x,y,\alpha} - d_aA_{x,y,\alpha} \\
        \frac{\partial P_{x,y}}{\partial t} &= qf_r(E_{x,y}P - d_pP_{x,y})
    \end{split}
    \label{IDC}
\end{align}

The eggs or immature populations at position $(x, y)$ in equation (\ref{IDC}) are grouped into sub-populations that have different $\alpha_i$ aggregation levels. New eggs are generated according to the stochastic process in Equation (\ref{povip}) and integrated into the existing population at $(x, y)$ with aggregation level $\alpha$, based on the number of adults with varying $\alpha$. The aggregation levels range from $\alpha = \alpha_n$ (minimum) to $\alpha = \alpha_x$ (maximum). The number of adults at $(x, y, \alpha)$ follows the process of the original compartmental model (\ref{compartmental}), and the number of natural enemies at $(x, y)$ is determined by one of the functional response models (Equations \ref{fr1}-\ref{fr3}). The term $f_r(E_{x,y,\alpha}P)$ represents the predator functional response, which can be one of the three functional response types presented in Equations \ref{fr1}-\ref{fr3}.

\paragraph{Boundary conditions}
To mitigate border effects, periodic boundary conditions were applied. These conditions assume that the values at the borders wrap around to the opposite side, as shown by the following equations:

\begin{equation}
    \begin{aligned}
        E(x, y) &= E(x + L_x, y) = E(x - L_x, y) \\
        E(x, y) &= E(x, y + L_y) = E(x, y - L_y) \\
        A(x, y) &= A(x + L_x, y) = A(x - L_x, y) \\
        A(x, y) &= A(x, y + L_y) = A(x, y - L_y) \\
        P(x, y) &= P(x + L_x, y) = P(x - L_x, y) \\
        P(x, y) &= P(x, y + L_y) = P(x, y - L_y)
    \end{aligned}
    \label{boundaryc}
\end{equation}

where $L_x$ and $L_y$ are the lengths of the matrix in the $x$ and $y$ directions, respectively.

\subsection{Model analysis}
The model analysis was conducted in three stages: \textit{a)} model dynamics, \textit{b)} allele distribution, and \textit{c)} spatial distribution of alleles.

\subsubsection{Numerical methods}

The system's properties were explored through numerical simulations of the spatio-temporal system. Descriptive statistics were calculated for state variable distributions, while time series analysis was carried out on the simulation output.

The system was solved using a hybrid approach combining the Finite Difference Method (FDM) \citep{hoffman2018numerical} for local reaction terms and spectral methods \citep{hesthaven2007spectral} via the Fast Fourier Transform (FFT) for diffusion. This FDM-FFT method effectively simulates reaction-diffusion processes, with FDM handling local reactions and FFT-based spectral methods managing diffusion.

Local reactions were computed pointwise at each grid location using finite differences, while diffusion was modelled by convolving with a diffusion kernel (e.g., Gaussian kernel) in Fourier space using FFT. This spectral approach improves computational efficiency and reduces numerical diffusion, especially on fine grids. The system was discretised into uniform grids in spatial, temporal, and genetic domains. The solution was calculated through the Euler scheme.

Preliminary simulations identified a representative combination of parameters where both natural enemies and prey persisted, using a type-I functional response and $\alpha = 1.0$ as the reference. The parameters for subsequent simulations are detailed in Table \ref{TP}.

\begin{table}[htpb]
\centering
    \small
    \begin{tabular}{l|cc}
    \toprule
    Parameter & Symbol & Value \\
    \midrule
    Population Model: & & \\
    New natural enemies per consumed prey/host & $q$ & 1.0 \\
    Hatching rate & $e$ & 1/2 \\
    Oviposition rate & $k$ & 3/4 \\
    Survival rate & $g$ & 1/2 \\
    Diffusion Coefficient & $D$ & 1/8 \\
    \midrule
    Functional response: & & \\
    Attack Rate FR I & $a$ & 1/4 \\
    Attack Rate FR II & $a$ & 1/4 \\
    Attack Rate FR III & $a$ & 1/4 \\
    Handling time & $h$ & 1.0 \\
    \midrule
    Egg model: & & \\
    $\alpha$ range ($\log$ scale) & $\alpha$ & {$-1,1$} \\
    Variability in $\alpha$ Heritability & $\sigma_h$ & 1/8 \\
    \midrule
    Simulation parameters: & & \\
    Matrix size & $x, y$ & $256, 256$ \\
    Time step & $dt$ & 1/32 \\
    Time length & $T$ & 1024 \\
    \bottomrule
    \end{tabular}
\caption{\label{TP} Parameters used for model simulations, along with their symbols and values.}
\end{table}

The $\alpha$ variable was defined as a vector of 17 values on a logarithmic scale, ranging from $\alpha = 0.1$ (near-equal egg masses) to $\alpha = 10.0$ (few large egg masses), covering $\alpha = \{0.1-10.0\}$. Simulations included four types of aggregation behaviour and three levels of natural enemies' functional responses, with $\alpha$ fixed at values 0.1, 1.0, and 10.0. The $\alpha$ parameter was varied over the range $\{0.1-10.0\}$, producing 12 combinations of functional responses and $\alpha$ levels. For each combination, random initial conditions were set, and the abundance of natural enemies ($P$), herbivores ($H$), and eggs ($E$) was recorded. Data were organised into five-dimensional matrices for herbivores $(x, y, t, \alpha, s)$ and four-dimensional matrices for natural enemies $(x, y, t, s)$.

Results were further used to calculate descriptive statistics and carry out time series analysis, exploring time dependencies and lags among state variables. The spatial domain ($x, y$) was represented by a $256 \times 256$ grid, with Equation (\ref{IDC}) iterated over 1024 time steps ($t = 1024$) through a time step of $1/64$. A hybrid approach combining Finite Difference Method (FDM) \citep{logan2022} was applied to reach numerical solutions. The diffusion term used a $3 \times 3$ convolution kernel in FFT space for faster processing. Integration in time variable was carried out using Euler's method \citep{butcher2016numerical}, selected for its speed and accuracy compared to the Runge-Kutta 4 method. Periodic boundary conditions were applied to mitigate border effects.

After preliminary tests to ensure numerical stability, we selected parameters shown in the final rows of Table \ref{TP}, using matrix sizes that minimise arithmetic errors. The model state was saved every 32 steps for analysis.

\paragraph{Initial conditions}

Initial conditions were randomly distributed using a gamma distribution, with parameters shown in Table \ref{TP}. The number of eggs and adult herbivores for each $\alpha$ aggregation level was uniformly distributed across the matrix, ensuring that observed variation in $\alpha$ was attributable to natural enemy responses, not initial spatial distributions.

\paragraph{Statistical analysis}

Time series data for the state variables were analysed to assess the relationship between allele abundance and total allele abundance. The cross-correlation function was computed for the last 512 time units using the \texttt{ccf} function from the \texttt{statsmodels} package \citep{seabold2010statsmodels}. Median and 95\% quantiles were considered to assess significant differences among functional response types or egg aggregation strategies.


\section{Results}

\noindent
Simulations showed distinct dynamics affected by combinations of egg aggregation behaviour and predator functional responses. The presentation of the results is organised as follows:

\begin{itemize}
\item[i)] Dynamics under an intermediate aggregation level ($\alpha = 1.0$) across three predator functional responses, representing the most common scenario in theoretical studies and serving as a reference for subsequent comparisons.
\item[ii)] Stability of the system across twelve combinations of three predator functional responses and four aggregation levels: extremely low ($\alpha = 0.1$), intermediate ($\alpha = 1.0$), extremely high ($\alpha = 10.0$), and variable, where aggregation behaviour evolves freely as described in Equation \ref{IDC}.
\item[iii)] Evolution of allele frequency for aggregation behaviour, identifying advantageous strategies in different scenarios and their impact on system dynamics.
\end{itemize}

\subsection{Model dynamics under intermediate aggregation levels}

\noindent
At intermediate aggregation levels with a type I functional response, the model showed classical Lotka-Volterra dynamics, characterised by periodic fluctuations of the of natural enemies ($P$) and prey/host ($H$) abundances. 
After the random initial distribution, the cycles maintained consistent periodicity, gradually decreasing in amplitude until stabilising at a mean abundance for all state variables (Fig. \ref{map_LVI}). The system stabilised after 1024 time units, with populations of $H$ and $P$ present in all simulations. The final spatial distribution resembled a random arrangement with uniform mean and variance.

Spatial and temporal patterns were significantly different in the case of type II functional responses. At the end of the simulations, the distribution of the abundance showed the wave pattern typical of spatially explicit Lotka-Volterra models, where the waves of $H$ were followed by waves of $P$, as observed in previous studies \citep{hastings1977spatial, hastings1978global, sherratt2001periodic, fu2015wave}. Oscillations were initially irregular, but a great portion of the simulations reached an equilibrium (Fig. \ref{map_LVII}). Type II responses also reported a small number of extinctions (Tables \ref{tab_res_H}, \ref{tab_res_E}, and \ref{tab_res_P}).

Spatial and temporal patterns were similar to type I models (Fig. \ref{map_LVIII}) in case of type III functional responses.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.0\textwidth]{Dynamics_Maps_01_00.png}
  \caption{\label{map_LVI} Dynamics of the compartmental model for herbivores with an egg predator (Equation \ref{IDC}) under a type-I functional response. Abundance maps: (A) Eggs ($E$), (B) Adult Herbivores ($H$), (C) Parasitoids/Predators ($P$); (D) Time series of mean abundances for all three populations.}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.0\textwidth]{Dynamics_Maps_02_00.png}
  \caption{\label{map_LVII} Dynamics of the compartmental model for herbivores with an egg predator (Equation \ref{IDC}) under a type-II functional response. Abundance maps: (A) Eggs ($E$), (B) Adult Herbivores ($H$), (C) Parasitoids/Predators ($P$); (D) Time series of mean abundances for all three populations.}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.0\textwidth]{Dynamics_Maps_03_00.png}
  \caption{\label{map_LVIII} Dynamics of the compartmental model for herbivores with an egg predator (Equation \ref{IDC}) under a type-III functional response. Abundance maps: (A) Eggs ($E$), (B) Adult Herbivores ($H$), (C) Parasitoids/Predators ($P$); (D) Time series of mean abundances for all three populations.}
\end{figure}


\clearpage
\subsection{Stability}

\noindent
Simulations for all combinations of functional responses and aggregation types (low, intermediate, high, and free-to-evolve) showed distinct dynamics. In case of type I functional responses simulations converged to a similar abundance level after dampened oscillations (Fig. \ref{map_LVOA}A-D). Lower aggregation levels ($\alpha = 0.1$) had minimal impact (Fig. \ref{map_LVOA}A, B), while higher $\alpha$ values ($\alpha = 10$) caused more irregular oscillations but resulted in lower equilibrium levels (Fig. \ref{map_LVOA}C). The free-to-evolve model produced dynamics similar to low aggregation levels but with slightly reduced oscillation amplitudes (Fig. \ref{map_LVOA}D). No extinctions occurred, indicating stability across scenarios (Tables \ref{tab_res_H}, \ref{tab_res_E}, and \ref{tab_res_P}).

In case of type II functional responses, instability at low aggregation levels was flanked by a higher extinction frequency. As aggregation increased, extinction rates decreased (Tables \ref{tab_res_H}, \ref{tab_res_E}, and \ref{tab_res_P}). At high aggregation levels or with free-to-evolve $\alpha$, no extinctions occurred, and irregular predator-prey cycles emerged, with increasing cycle amplitudes (Fig. \ref{map_LVOA}F-H).

In case of type III functional responses, the system remained stable without extinctions (Tables \ref{tab_res_H}, \ref{tab_res_E}, and \ref{tab_res_P}). At intermediate-to-low aggregation levels, dampened oscillations led to equilibrium (Fig. \ref{map_LVOA}I, J), similar to type-I dynamics. Higher $\alpha$ levels produced irregular oscillations with increased amplitude, while free-to-evolve $\alpha$ led to even higher amplitude aggregations over time, without extinctions.

\begin{table}[htpb]
    \centering
    \small
    \begin{tabular}{ll|cccc}
     FR type & Statistic & $\alpha = 0.1$ & $\alpha = 1.0$ & $\alpha = 10.0$ & $\alpha = \{0.1-10.0\}$\\
     \hline
     I       & EF        & 0.00           & 0.00           & 0.00            & 0.00 \\
             & M         & 0.60           & 0.60           & 0.42            & 0.60\\
             & Q         & (0.60, 0.60)    & (0.60, 0.60)    & (0.42, 0.43)     & (0.60, 0.60)\\
     II      & EF        & 0.75           & 0.08           & 0.00            & 0.00 \\
             & M         & $-\infty$      & 0.67           & 1.54            & 3.68\\
             & Q         & ($-\infty$, 0.62) & ($-\infty$, 0.76) & (1.33, 2.36)   & (3.07, 4.83)\\
     III     & EF        & 0.00           & 0.00           & 0.00            & 0.00\\
             & M         & 0.53           & 0.53           & 1.21            & 2.63\\
             & Q         & (0.53, 0.53)    & (0.53, 0.53)    & (1.06, 1.82)    & (1.63, 4.27)\\
     \hline
    \end{tabular}
    \caption{Extinction frequency (EF), median (M), and the 95\% quantiles (Q) of the abundances of the global population sizes of Herbivores at the end of the simulations. The abundances are expressed in base-10 logarithms.}
    \label{tab_res_H}
\end{table}

\begin{table}[htpb]
    \centering
    \small
    \begin{tabular}{ll|cccc}
     FR type & Statistic & $\alpha = 0.1$ & $\alpha = 1.0$ & $\alpha = 10.0$ & $\alpha = \{0.1-10.0\}$\\
     \hline
     I       & EF        & 0.00           & 0.00           & 0.00            & 0.00 \\
             & M         & 0.60           & 0.60           & 0.42            & 0.60\\
             & Q         & (0.60, 0.60)    & (0.60, 0.60)    & (0.42, 0.43)     & (0.60, 0.60)\\
     II      & EF        & 0.75           & 0.08           & 0.00            & 0.00 \\
             & M         & $-\infty$      & 0.67           & 1.54            & 3.63\\
             & Q         & ($-\infty$, 0.62) & ($-\infty$, 0.76) & (1.32, 2.4)   & (2.09, 4.76)\\
     III     & EF        & 0.00           & 0.00           & 0.00            & 0.00\\
             & M         & 0.54           & 0.54           & 1.20            & 2.57\\
             & Q         & (0.54, 0.54)    & (0.54, 0.54)    & (1.05, 1.76)    & (1.58, 4.21)\\
     \hline
    \end{tabular}
    \caption{Extinction frequency (EF), median (M), and the 95\% quantiles (Q) of the abundances of the global population sizes of Eggs at the end of the simulations. The abundances are expressed in base-10 logarithms.}
    \label{tab_res_E}
\end{table}

\begin{table}[htpb]
    \centering
    \small
    \begin{tabular}{ll|cccc}
     FR type & Statistic & $\alpha = 0.1$ & $\alpha = 1.0$ & $\alpha = 10.0$ & $\alpha = \{0.1-10.0\}$\\
     \hline
     I       & EF        & 0.00           & 0.00           & 0.00            & 0.00 \\
             & M         & 0.60           & 0.60           & 0.39            & 0.59\\
             & Q         & (0.60, 0.60)    & (0.60, 0.60)    & (0.38, 0.40)     & (0.59, 0.59)\\
     II      & EF        & 0.75           & 0.08           & 0.00            & 0.00 \\
             & M         & $-\infty$      & 0.67           & 1.51            & 5.18\\
             & Q         & ($-\infty$, 0.59) & ($-\infty$, 0.75) & (1.30, 2.23)   & (3.56, 6.25)\\
     III     & EF        & 0.00           & 0.00           & 0.00            & 0.00\\
             & M         & 0.53           & 0.53           & 1.21            & 3.32\\
             & Q         & (0.53, 0.53)    & (0.53, 0.53)    & (1.05, 2.01)    & (1.60, 5.95)\\
     \hline
    \end{tabular}
    \caption{Extinction frequency (EF), median (M), and the 95\% quantiles (Q) of the abundances of the global population sizes of Natural Enemies at the end of the simulations. The abundances are expressed in base-10 logarithms.}
    \label{tab_res_P}
\end{table}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.0\textwidth]{LVOA_fr_3_a_02.png}
  \caption{\label{map_LVOA} Time series for the egg aggregation model with three aggregation levels and the free-to-evolve model, incorporating three functional response types. Each row represents a different functional response, with the first row showing type-I, the second row type-II, and the third row type-III. The columns represent different aggregation levels ($\alpha$), with Natural Enemies abundance in green, Eggs abundance in blue, and Adults abundance in red.}
\end{figure}

\clearpage

\subsection{Allele Frequency}
\noindent
The allele frequencies over time varied across the three types of predator functional responses. In case of type I functional responses, characterised by stable dampened oscillations, allele frequency variation was smooth. Alleles linked to low aggregation gradually replaced the alleles associated with high aggregation; alleles with $\alpha > 3.16$ extinguished after 200 simulation time units in all scenarios (Fig. \ref{map_LV2DT}A, D). Consequently frequencies of high $\alpha$ alleles approached zero in the end of the simulation (Fig. \ref{map_LV2DH}A).

The cross-correlation function (Fig. \ref{ccf_a}A) showed that the $\alpha = 0.1$ allele lagged slightly behind total herbivore abundance, with minimal difference in lag. A sinusoidal pattern showed that this allele depended on herbivore abundance oscillations. The allele $\alpha = 1.0$ showed a smaller lag, while for $\alpha = 10.0$ the population extinguished early on, with only sporadic appearances due to random mutations (Equation \ref{hovip}).

In case of type II predator functional responses, the pattern reversed. At the beginning, alleles with high $\alpha$ values decreased in abundance, as in type I responses, until population bottlenecks occurred around simulation time 32. Afterwards, low $\alpha$ alleles largely extinguished (Fig. \ref{map_LV2DT}B, E), leaving only alleles with $\alpha > 1$ in low frequencies until the end (Fig. \ref{map_LV2DH}B). Allele frequencies became highly variable following bottlenecks, driven by the competitive advantage of lower $\alpha$ alleles outside these bottlenecks. High-alpha alleles showed declining trends, punctuated by sudden increases caused by local bottlenecks from predator-prey waves.

Under type II responses, alleles with $\alpha = 0.1$ and $\alpha = 1.0$ slightly preceded herbivore abundance, while $\alpha = 10.0$ trailed it (Fig. \ref{ccf_a}B). This suggests that after each herbivore bottleneck, intermediate and low aggregation alleles extinguished, recovering through mutation (Equation \ref{hovip}).

In case of type III functional responses, allele distribution dynamics resembled those of type II, but with smoother transitions. High $\alpha$ alleles gradually increased in abundance, while low $\alpha$ alleles never completely disappeared (Fig. \ref{map_LV2DT}C, F). Consequently, low-alpha alleles remained at low levels throughout simulations, never reaching zero (Fig. \ref{map_LV2DH}C). The relative abundance of intermediate to high aggregation alleles ($\alpha = 1.0$ and $\alpha = 10.0$) trailed herbivore abundance, while the $\alpha = 0.1$ allele did not show a clear pattern, although its cross-correlation function showed a peak preceding herbivore abundance (Fig. \ref{ccf_a}C). This indicates that the $\alpha = 0.1$ allele extinguished after each bottleneck, only to recover slightly through mutation, remaining at very low levels and resulting in noisy, non-significant cross-correlation patterns.

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.0\textwidth]{FR_time_alleles.png}\\
  %\includegraphics[width=1.0\textwidth]{FR_time_alleles_img.png}
  \caption{\label{map_LV2DT} Time series of allele frequencies for all functional response types. Panels A, B, and C show the median and 95\% quantiles of allele frequencies for the lowest ($\alpha = 0.1$), intermediate ($\alpha = 1.0$), and highest aggregation ($\alpha = 10.0$), respectively, for functional responses I, II, and III. The lower row displays the mean allele frequencies across all simulations.}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.1\textwidth]{Alleles_hist.png}
  \caption{\label{map_LV2DH} Final histogram of allele frequency distributions across all simulations for the herbivorous populations. The lines represent the percentage of simulations where each allele had the frequency indicated by the x-axis. Alleles are represented by different colours: purple for low aggregation, red for high aggregation, and green for intermediate values.}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[width=1.1\textwidth]{Alleles_ccf.png}
  \caption{\label{ccf_a} Cross-correlation function between the relative frequency of alleles associated with different aggregation levels: lowest ($\alpha = 0.1$, in red), intermediate ($\alpha = 1.0$, in green), and highest ($\alpha = 10.0$, in blue), and total herbivore abundance. Positive lag values indicate when herbivore abundance lags behind allele abundance, while negative lag values indicate the opposite. Solid lines represent the median, with shaded areas representing the 95\% quantiles. The solid circle marks the peak of each cross-correlation function.}
\end{figure}
\clearpage

\section{Discussion}

\noindent
This study explored, from a modelling perspective, egg aggregation as an evolutionarily stable strategy for avoiding specific natural enemies and the effects of reduced enemy efficiency at high host densities. These aspects are critical in ecological and pest control contexts, especially in managing outbreaks of species that are difficult to contain and cause damage in cultivated fields or natural ecosystems, such as forests. Biological control, using natural enemies instead of chemical pesticides, is increasingly important as pesticide use faces growing restrictions \citep{DiSora2024}, or is prohibited at all \citep{BonoRossello2023}. Previous research has focused on host-prey and parasitoid-predator interactions to predict the optimal timing for releasing natural enemies \citep{van2010quality, singh2016aphids, ong2015coupling}. However, the evolutionary and behavioural responses of prey/hosts must also be considered, as these mechanisms are crucial for species survival and long-term persistence.

This study highlights the substantial role of spatial egg aggregation in population dynamics and system stability, although the optimal oviposition strategy depends on the nature of the biological controllers. Egg aggregation offers advantages, but it is not universally adopted by the species; its effectiveness is sometimes constrained by substrate or environmental factors \citep{rohlfs2004spatial}. Furthermore, some species show contrasting strategies: early larval instars aggregate to evade natural enemies, whereas later instars disperse in response to the same threats \citep{cornell1987developmental}.

At high population densities, aggregation patterns may reverse or become overdispersed \citep{battisti2023aggregative}, suggesting that a species’ strategy can shift depending on the environmental context. As noted by \cite{hassell1991spatial}, local movement in spatially explicit Lotka-Volterra models can enhance the persistence of unstable populations. Our results support this piece of knowledge, with some scenarios that show spatial patterns akin to random white noise.

Factors such as environmental noise colour \citep{bruzzone2023effect} or delayed egg emergence \citep{corley2009delayed} can trigger insect outbreaks. Some researchers \citep{keeling2000reinterpreting} suggest that delays, functional responses, and spatial heterogeneity are similar phenomena in population dynamics, leading to comparable outbreak patterns and spatial distributions.

Under type I predator functional responses, the system remained stable, with equilibrium values independent of egg aggregation. This evidence suggests that aggregation has little impact on the effectiveness of natural enemies. In contrast, type II responses showed that egg aggregation supported herbivore populations to surviving outbreaks, however, they paradoxically destabilised the system leading to boom-and-bust cycles. Even highly aggregated strategies could not prevent extinction under type II dynamics. Type III functional responses showed a greater stability than type II, with low prey densities rendering natural enemies ineffective, thus preventing extinction. At high densities, a refuge effect emerged, akin to type II dynamics, where low $\alpha$ levels resembled type I dynamics and high $\alpha$ levels resembled type II dynamics.

In type I systems, predation or parasitism is independent on prey/host density, meaning that egg aggregation provides no advantage and results in lower $\alpha$ values. However, in type II and III systems, functional responses reach a plateau at high densities, creating inverse density dependence and a refuge effect. In these systems, aggregating oviposition provides a competitive advantage, supporting \cite{chesson1986aggregation}'s aggregation-of-risk hypothesis.

The results of this study support the notion that minimal aggregation spreads risk across space, while increased aggregation concentrates risk in specific areas. In low aggregation, natural enemies struggle to find isolated eggs, incurring in higher search costs. In contrast, high aggregation leads to risk dilution, as predators are less efficient at handling large quantities of prey, allowing some eggs to evade predation \citep{desurmont2011aggregative, chesson1986aggregation}.

This study is a theoretical exploration within a uniform area with unlimited resources for a single herbivore species targeted by one natural enemy. While the setup may seem focused, the results align with existing research on outbreaking and non-outbreaking species. Studies such as \citet{cappuccino1995population, stamp1980egg, hunter1991traits} indicate that while many species spread oviposition to evade generalist predators, specific natural enemies exert significant selective pressure.

Future research should incorporate more complex models, including heterogeneous landscapes, limited resources, refuges, or field data feedback \citep{van2017spatially, bruzzone2023effect, BonoRossello2023, welch2014temporal}. Experimental setups in controlled conditions could help validate the hypotheses generated by these models.

Our findings suggest a reinforcing cycle between evading specific natural enemies and the frequency of outbreaks. Aggregative oviposition helps avoid density-dependent predation or parasitism, favouring aggregative alleles and increasing outbreak likelihood. This feedback loop could present significant challenges in pest management, as aggregative behaviour reduces predator efficiency, allowing herbivore populations to grow and overwhelm plant defences, thus exacerbating outbreaks \citep{wetzel2023variability}.

Furthermore, systems with type II functional responses may experience even more frequent outbreaks than those with type I responses, as the reinforcing cycle is stronger in such systems.

\section{Conclusions}
\noindent
In conclusion, egg aggregation is an important strategy for evading specific natural enemies under type II and type III systems but has little effect under type I systems. The resulting instability in type II and type III systems enhances the selective pressure for aggregative oviposition, further amplifying outbreak frequency in a positive feedback loop.

\section*{Acknowledgements}
\noindent
We sincerely thank Carlos Bernstein for insightful discussions on the theoretical framework of an early version of this model in 2005. While his suggestions were not directly incorporated into the final study, they significantly shaped the conceptual development and interpretation of our results. We also extend our gratitude to Guillermo Logarzo for his invaluable guidance in understanding the functional response mechanisms of predators and parasitoids, which greatly contributed to refining our model.

\section*{Funding}
\noindent
Luca Rossini is supported by the European Commission under the Marie Sklodowska Curie Actions Postdoctoral Fellowship (MSCA-PF-2022) project ``PestFinder'' (Grant No. 101102281). José Crespo is funded by the National Scientific and Technical Research Council of Argentina through the PIP 1939 grant.

\section*{Conflict of Interest}
\noindent
The authors declare no conflict of interest.

\section*{Ethics Approval}
\noindent 
Ethical approval was not required for this study, as it consists solely of simulations of theoretical populations without reference to any specific geographic location or ecosystem.

\section*{Consent to Participate}
\noindent
Not applicable, as this is a theoretical study.

\section*{Consent for Publication}
\noindent
Not applicable, as this is a theoretical study.

\section*{Availability of Data and Material}
\noindent
Not applicable, as this study is based entirely on simulations.

\section*{Code Availability}
\noindent
The code used for simulations, statistical analyses, and visualizations is available at:
\texttt{https://github.com/okktawio/Egg-Aggregation-Outbreaks}.

\section*{Author Contributions}
\noindent
O.A. Bruzzone (OAB) conceptualized the study, developed the model, and conducted the analysis. L. Rossini (LR) reviewed the mathematical framework, results, and main text. J. Crespo (JC) contributed to the conceptualization and reviewed the mathematical aspects. M.B. Aguirre (MBA) conducted bibliographic research, revised the manuscript, and assisted in organizing the results, tables, and figures. All authors contributed to interpreting the results, writing the manuscript, and reviewing the final version.

\section*{Declaration of Generative AI and AI-Assisted Technologies in the Writing Process}
\noindent
During manuscript preparation, the authors used OpenAI's ChatGPT 3.5 and Grammarly to refine wording and grammar. After using these tools, the authors reviewed and edited the content as needed and take full responsibility for the final version.

\bibliographystyle{apalike} 
\bibliography{egg_aggregation}

\end{document}
