A predator-prey model with age-structured role reversal.

A predator-prey model with age-structured role reversal - PMC Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice J Math Biol . 2026 Apr 17;92(5):73. doi: 10.1007/s00285-026-02402-5 Search in PMC Search in PubMed View in NLM Catalog Add to search A predator-prey model with age-structured role reversal Luis Carlos Suarez Luis Carlos Suarez 1 Department of Mathematics, University of Maryland, 4176 Campus Drive, College Park, 20742 MD USA Find articles by Luis Carlos Suarez 1, # , Maria K Cameron Maria K Cameron 1 Department of Mathematics, University of Maryland, 4176 Campus Drive, College Park, 20742 MD USA Find articles by Maria K Cameron 1, ✉, # , William F Fagan William F Fagan 2 Department of Biology, University of Maryland, 4094 Campus Drive, College Park, 20742 MD USA Find articles by William F Fagan 2 , Doron Levy Doron Levy 1 Department of Mathematics, University of Maryland, 4176 Campus Drive, College Park, 20742 MD USA Find articles by Doron Levy 1 Author information Article notes Copyright and License information 1 Department of Mathematics, University of Maryland, 4176 Campus Drive, College Park, 20742 MD USA 2 Department of Biology, University of Maryland, 4094 Campus Drive, College Park, 20742 MD USA ✉ Corresponding author. # Contributed equally. Received 2025 Feb 26; Revised 2026 Mar 24; Accepted 2026 Apr 8; Issue date 2026. © The Author(s) 2026 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ . PMC Copyright notice PMCID: PMC13090279  PMID: 41995875 Abstract We propose a predator-prey model with an age-structured predator population that exhibits a functional role reversal. The structure of the predator population in our model embodies the ecological concept of an "ontogenetic niche shift," in which a species’ functional role changes as it grows. This structure adds complexity to our model but increases its biological relevance. The time evolution of the age-structured predator population is motivated by the Kermack-McKendrick Renewal Equation (KMRE). Unlike KMRE, the predator population’s birth and death rate functions depend on the prey population’s size. We establish the existence, uniqueness, and positivity of the solutions to the proposed model’s initial value problem. The dynamical properties of the proposed model are investigated via Latin Hypercube Sampling in the 15-dimensional space of its parameters. Our Linear Discriminant Analysis suggests that the most influential parameters are the maturation age of the predator and the rate of consumption of juvenile predators by the prey. We carry out a detailed study of the long-term behavior of the proposed model as a function of these two parameters. In addition, we reduce the proposed age-structured model to ordinary and delayed differential equation (ODE and DDE) models. The comparison of the long-term behavior of the ODE, DDE, and the age-structured models with matching parameter settings shows that the age structure promotes the instability of the Coexistence Equilibrium and the emergence of the Coexistence Periodic Attractor. Keywords: Age-structured, Role reversal, Kermack-McKendrick renewal equation, Latin hypercube sampling, Linear discriminant analysis, Phase diagrams Introduction An overview People have been thinking about trophic interactions between various species since ancient times. The concept of a food chain was introduced by a 10th-century Arab philosopher, al-Jahiz (Agutter and Wheatley 2008 ). The birth of modern ecology, recognizing complex interactions between species and equipped with mathematical modeling, can be attributed to the 1920s, when Elton ( 1927 ) upgraded food chains to food webs and Lotka ( 1925 ) and Volterra ( 1926 ) proposed their predator-prey models. While food webs can sometimes correctly capture the big picture for ecological systems, they typically oversimplify complex inter- and intraspecific interactions by ignoring an aspect of crucial importance: the variance in body size among conspecific individuals. Such variance is often associated with age differences, but other factors, such as differential foraging success, can also lead to variance in body size. Werner and Gilliam ( 1984 ) argued that changes in the body size of animals have a big influence on ecological processes. For example, juvenile predators may have very different diets from conspecific adults, and they may be eaten by other species that also compete with the prey of conspecific adults for resources. Such interactions, which biologists term “ontogenetic niche shifts," mean that a species’ functional role changes as it grows. Such shifts may negatively affect the survival and recruitment of predators, creating the so-called juvenile predator bottleneck. Further complexity of ecosystems is caused by role-reversal and cannibalism (Polis 1988 ). Polis ( 1991 ) pointed out four problems commonly present in catalogs of empirical food webs, making them inadequate for analysis: aggregation and undercounting of species, inadequate dietary information, lack of attention to age structure and ontogenetic niche shifts, and neglect of ecological feeding loops that occur due to role reversal and cannibalism. Polis gave several examples of ontogenetic reversals in the ecosystem of the Coachella Valley in Riverside County, California. For instance, gopher snakes eat eggs and young burrowing owls, while adult burrowing owls eat young gopher snakes. In fact, biological instances of role reversal and cannibalism are quite common in nature (Polis 1988 , 1991 ). For example, in the small deciduous forest on the campus of Hokkaido University, Japan, spider mites kill the larvae of one of their key predators (Saitā 1986 ). In another system, cod, the top predator in the Baltic Sea, prey on herring and sprat that consume the eggs of cod (Köster and Möllmann 2000 ), and juvenile cod between ages 0 and 2 may also lose from about a third to nearly half of its population due to cannibalism (Neuenfeldt and Köster 2000 ). In some systems, these interactions may be sublethal, such as the size-dependent aggression towards siblings that takes place among Dendrobates tinctorius tadpoles (Fouilloux et al. 2022 ). In the review paper of 2015, Nakazawa ( 2015 ) gave examples of three types of ontogenetic niche shifts—in diet, interaction type, and habitat—and emphasized the necessity of accounting for ontogenetic perspectives for understanding mechanisms underlying biodiversity and ecosystem functioning. Mathematical modeling of ecosystems Numerous mathematical models of ecological systems, with various degrees of complexity, were designed to capture and analyze particular aspects of the systems’ dynamics. Models with age-structured predator population based on the McKendrick equation (Mckendrick and Pai 1912 ) were developed and studied by Cushing and Saleem ( 1982 ) and  Cushing ( 1984 , 1986 ). Size-structured models featuring cannibalism based on the continuity equation were proposed by Cushing ( 1992 ) and Fagan and Odell ( 1996 ). An ODE-based predator-prey model with role reversal developed by Lehtinen ( 2021 ) demonstrated the feasibility of many ecological scenarios. The sudden ecological shifts from coexistence to predator extinction may occur due to various catastrophic bifurcations, such as saddle-node, homoclinic, and subcritical Hopf. Li et al. ( 2022 ) proposed a simple and enlightening generic ODE-based prey-predator model with role reversal. We use their model as a building block in this work. Delayed differential equations (DDEs) offer an alternative and, perhaps, simpler way to account for the age structure than the McKendrick equation (Mohr et al. 2014 ). Recently, Mishra et al. ( 2024 ) introduced a DDE-based model with role reversal. They showed that the maturation time of the predator and the prey handling time are the key factors in the dynamics of this model. These models of complex interactions between predator and prey can be summarized as follows. The dynamics depend on model type (ODE, DDE, PDE, etc.), model settings, and model parameters and can be complex. The predictions of various models are hard to compare to each other as the authors focus on different aspects of ecological processes. The goal and a summary of main results Here, our goal is to understand how age structure of the predator population and role reversal jointly affect the dynamics of a coupled predator-prey system. The first objective is to develop and investigate a predator-prey model with an age-structured predator population and role reversal. The second objective is to examine how the model type affects the long-term dynamics. Thus, we develop a mixed ordinary and partial differential equation-based model featuring an age-structured predator population and role reversal. The time evolution of the age density of the predator population is modeled by a Kermack-McKendrick-type Renewal Equation with birth- and death-rate functions depending on the prey population size. Note that the term population size refers to the density or numbers of prey or predators rather than their physical size. The prey population, in turn, depends on the juvenile and adult predator population sizes. This makes the predator PDE analytically unsolvable in contrast to that in Mishra et al. ( 2024 ). The dynamics of the prey population are governed by an ODE mimicking the prey equation in Li et al. ( 2022 ). We analyze the proposed model using analytical and numerical tools. We prove the existence, uniqueness, and positivity of the solution to the initial value problem for the proposed model, and the convergence of the numerical solution to it. In addition, we obtain systems of ODEs and, separately, DDEs that approximate the dynamics of the age-structured predator population. The proposed model involves 15 parameters. We use the Latin Hypercube Sampling Method (Audze and Eglãjs 1977 ; McKay et al. 1979 ; Iman et al. 1981 ) to examine the long-term behavior of the model over the parameter space. The results show that the system settles to one of three attractors, which we term the Equilibrial Coexistence Attractor, the Periodic Coexistence Attractor, and the Predator-Free Attractor, in 22%, 19%, and 55% of cases, respectively. In the remaining approximately 4% of cases, the prey and predator population sizes blow up. This blow-up is possible due to unbounded birth rates of prey and predators and a complex interplay between their population sizes involving time delay effects. We conduct a Linear Discriminant Analysis (LDA) (a.k.a. Multiple Discriminant Analysis (MDA)) (Duda et al. 2001 ) and conclude that the two most influential parameters on the qualitative long-term behavior of the system are the maturation age τ ∗ and the consumption rate of juvenile predators by prey g . Subsequently, we fix most parameters at some default values and conduct a thorough investigation of the behavior of the proposed model on the two most pivotal parameters, τ ∗ and g . We plot phase diagrams in ( τ ∗ , g ) -space, denoting regions of the various attractor types, and a collection of bifurcation diagrams. We do so for two settings of the model characterized by a sharp or a gradual transition from the juvenile predator to adult. The regions of various attractor types in these settings are somewhat different. We present a collection of bar plots of the age density of the predator at Equilibrial Coexistence Attractors and at a Periodic Coexistence Attractor. The age density is a monotonically decreasing and fast-decaying function for the Equilibrial Coexistence Attractors. It exhibits traveling waves with decaying amplitude during the period of the Periodic Coexistence Attractor. In all cases, the predator population size decays to negligibly small values before reaching the chosen age cap. Finally, we derive ODE- and DDE-based models from our age-structured model, aiming to make the equations of all three models consistent. We compare the phase diagrams of these three models in the ( τ ∗ , g ) -space. Evidently, the age structure of the predator population promotes oscillatory behavior. To be precise, Periodic Coexistence Attractors exist in all three models, but the sizes of the regions in the parameter space ( τ ∗ , g ) where the Periodic Coexistence Attractor is observed have drastically different sizes across models: the largest for the age-structured model, the smallest for the ODE model, and of an intermediate size for the DDE model. Our codes for simulating the age-structured model and the corresponding ODE and DDE models and plotting phase and bifurcation diagrams are available on GitHub (Suarez and Cameron 2025 ). The remainder of the paper is organized as follows. The age-structured model is developed in Section 2. The conditions for existence, uniqueness, and positivity are stated in Section 3 and proven in Appendix A . The ODE and DDE models are derived in Sections 4 and 5 respectively. Numerical algorithms are described in Section 6 . The results of numerical studies of the age-structured model and the corresponding ODE and DDE models are presented in Section 7 . Our findings are discussed in Section 8 , and conclusions are drawn in Section 9 . Model development In this section, we develop a mixed ODE- and PDE-based predator-prey model with age-structured predator population and role reversal. A building block: Li, Liu, and Wei’s ODE Li et al. ( 2022 ) proposed the following ODE model for a generic predator-prey system with role reversal: x ′ = x ( r - a x + s y 1 - b y 2 ) , y 1 ′ = k x y 2 - y 1 ( g x + D + m 1 ) , y 2 ′ = D y 1 - m 2 y 2 . 1 In this model, the population of predators is subdivided into juvenile and adult predators. x , y 1 , and y 2 represent the population sizes of prey and the juvenile and adult predator respectively, while a , s , b , k , g , D , m 1 and m 2 are parameters whose biological meaning and values are listed in Table 1 – we borrowed the selected values from Li et al. ( 2022 ) except for that of b that we increased from 0.4 in (Li et al. 2022 ) to 0.8 to obtain richer dynamics. Juvenile predators are assumed to lack hunting skills and depend on adult predators for food. However, due to their assumed smaller size, juvenile predators may be eaten by the prey species, which also obtains resources from the environment and is, in turn, eaten by adult predators. Table 1. Parameters of the proposed model. The values of r , a , k , and s are borrowed from Li et al. ( 2022 ) Parameter Biological meaning Selected Value Range τ ∗ Maturation age of predator 1 [0,2] g Consumption rate of juvenile predator by prey 0.2 [0,1] ν Smoothness of the indicator function 1 and 100 [1,100] r Intrinsic growth rate of prey 0.4 [0.1,0.6] a Intraspecific competition rate of prey 0.01 [0.005,0.05] k Reproduction rate of predator 0.3 [0.1,1] b Consumption rate of prey by adult predator 0.8 [0.1,1] s Growth rate of prey due to eating juvenile predators 0.2 [0, 1, 1] ζ Birth rate, exponential factor 10 [5,20] μ M Prefactor for the hunger death rate of predator 1 [0.5,5] ρ Exponential factor of the hunger death rate of predator 5 [3,7] d p Prefactor of the base death rate of predator 0.4 [0.1,1] b p Prefactor of the base birth rate of predator 0.05 [0.03,0.1] b ep Exponential factor of the base birth rate of predator 0.1 [0.05,0.15] d ep Exponential factor of the base death rate of predator 0.1 [0.05,0.15] Open in a new tab The trivial extinction equilibrium E 0 = ( 0 , 0 , 0 ) is always a saddle point with a single positive eigenvalue r corresponding to the eigenvector (1, 0, 0). The other equilibria of ( 1 ), the Predator-Free Equilibrium and the Coexistence Equilibrium E ⋆ = ( x ⋆ , y 1 ⋆ , y 2 ⋆ ) , are asymptotically stable or unstable depending on parameter values. When E ⋆ = ( x ⋆ , y 1 ⋆ , y 2 ⋆ ) is unstable, there exists the Periodic Coexistence Attractor around it Li et al. ( 2022 ). The important conditions for the existence of attractors where both species maintain population sizes bounded away from zero are s ≤ g and k ≤ b (Li et al. 2022 ). These conditions mean that the eaten population is damaged more than the eating population benefits from consuming it. While model ( 1 ) captures that juvenile predators do not reproduce and lack hunting skills, it is still highly idealized. Most importantly, it ignores the time delay between the birth of juvenile predators and their maturation. Instead, juvenile predators become adults at the rate D . The Kermack-McKendrick Renewal Equation The dynamics of an age-structured population are described by the classical Kermack-McKendrick Renewal Equation (KMRE) (Perthame 2006 , Section 3.1) u t ( t , τ ) + u τ ( t , τ ) = - μ ( τ ) u ( t , τ ) , u ( t , 0 ) = ∫ 0 ∞ B ( τ ) u ( t , τ ) d τ , u ( 0 , τ ) = u 0 ( τ ) . 2 Here, the independent variables t and τ represent time and age, respectively; u ( t , τ ) is the age density, and μ ( τ ) and B ( τ ) are the death and birth rate functions, respectively. The age of the generation born at time t 0 is τ = t - t 0 . KMRE says that the age density of any generation decreases at the age-dependent rate μ ( τ ) . The initial age density of the predator population, u 0 ( τ ) , is bounded and integrable. Furthermore, the birth and death rate functions are bounded and satisfy the condition 1 < ∫ 0 ∞ B ( τ ) e - M ( τ ) d τ < ∞ , where M ( τ ) : = ∫ 0 τ μ ( s ) d s . 3 These assumptions are enough to guarantee the existence and uniqueness of the solutions (Perthame 2006 ). Furthermore, equation ( 2 ) admits an implicit analytical solution under these assumptions, obtained via the method of characteristics: u ( t , t + τ 0 ) = u 0 ( τ 0 ) e - ∫ 0 t μ ( t ′ + τ 0 ) d t ′ , τ ≡ t + τ 0 ≥ t , u ( t 0 , t 0 + τ ) = ∫ 0 ∞ B ( τ ′ ) u ( t 0 , τ ′ ) d τ ′ e - ∫ 0 τ μ ( τ ′ ) d τ ′ , t ≡ t 0 + τ ≥ τ , 4 where τ 0 is the age at time zero and t 0 is the birth time. An age-structured predator-prey PDE-based model We will model the predator population by a renewal equation similar to KMRE ( 2 ). The difference between KMRE and our renewal equation governing the predator population dynamics with the role reversal is two-fold. First, the birth and death rate functions of the predator depend on the prey population size, x , governed by an ODE. Second, the predator population is divided into two subpopulations of juvenile and adult predators, where the prey can eat the juvenile predators and are, in turn, eaten by the adult predators. An important parameter of our model is the maturation age τ ∗ > 0 . It is used to split the total predator population, y ( t ), to juvenile, y 1 ( t ) , and adult, y 2 ( t ) , populations and to define predator’s birth and death rate functions, B ( x , τ ) and μ ( x , τ ) . The juvenile, adult, and total predator populations are given, respectively, by y 1 ( t ) = ∫ 0 τ ∗ u ( t , τ ) d τ , y 2 ( t ) = ∫ τ ∗ ∞ u ( t , τ ) d τ , y ( t ) = y 1 ( t ) + y 2 ( t ) . The ODE for the prey population, x , is borrowed from Li et al. ( 2022 ). The predator birth and death rate functions, B and μ , are chosen to be of the form B ( x , τ ) = k x φ { τ ≥ τ ∗ } ( τ ) + B ~ ( τ ) ( 1 - e - ζ x ) , μ ( x , τ ) = g x φ { τ < τ ∗ } ( τ ) + μ B ( τ ) + μ M e - ρ x . Here, ζ and ρ are positive constants. The functions φ { · } are smooth approximations of the characteristic functions of [ τ ∗ , ∞ ) and ( - ∞ , τ ∗ ] defined, respectively, by φ { τ ≥ τ ∗ } ( τ ) = 1 1 + e - ν ( τ - τ ∗ ) , φ { τ < τ ∗ } ( τ ) = 1 1 + e - ν ( τ ∗ - τ ) where ν is a positive constant. As ν → ∞ , the functions φ { · } converge pointwise to the indicator functions of the corresponding intervals. As ν → 0 , φ { · } converges pointwise to for all τ . We chose φ { · } to be smooth to account for the gradual maturation process that is typical of many species. In nature, the transition from juveniles to adults may span a period from a few minutes to many years, depending on the species. The function B ~ ( τ ) is the base birth rate supported on [ τ ∗ , ∞ ) so that juvenile predators are unable to reproduce. It is nonnegative, bounded, and integrable. Importantly, the birth rate function B ( x , τ ) is zero at any age τ if the prey population size x is zero, and it is positive at any positive x . The choice of the death rate function μ ( x , τ ) is motivated by the expected behavior of the predator. It is the sum of three terms describing three reasons for predators to die. The juvenile predator death rate, g x φ { τ < τ ∗ } , is due to being eaten by prey. The base predator death rate, μ B ( τ ) , accounts for deaths for age-related reasons. The term μ M e - ρ x describes the death rate from hunger. In summary, the proposed model with the age-structured predator population and role reversal is given by The initial size of the prey population x 0 is positive. The initial age density of the predator, u 0 ( τ ) , is nonnegative and has a compact support. Existence, uniqueness, and positivity In this section, we will establish the existence, uniqueness, and positivity of the solution to the system (5a)–(5g). These properties are naturally expected from a model describing population dynamics. However, the task of developing a model with these properties is nontrivial. Before settling on the model (5a)–(5g), we designed a different model motivated by the ODE model ( 1 ) from Li et al. ( 2022 ) but abandoned it because it led to a U-shaped age density of the predator. Furthermore, proving these properties for the system (5a)–(5g) turned out to be harder than we expected. Similar properties for ( 1 ) directly follow from Nagumo-Brezis’ Theorem (Theorem 1.2, page 26 in Pavel and Motreanu ( 1999 )) about flow invariance. Unfortunately, this theorem is not applicable to the system (5a)–(5g) because the flow is not of the form z ′ = F ( t , z ) due to the renewal of the characteristics. The proof of existence, uniqueness, and positivity of the solution to the Kermack-McKendrick Renewal Equation in (Iannelli 1995 )(Chapter VII) exploits the fact that the solution can be written out analytically in the form of integrals. This is impossible to do for our model (5a)–(5g) because of the interdependence of x and u as the predator’s birth and death rates, B ( x , t ) and μ ( x , t ) , depend on x , and x depends on u . As a result, our proof of the basic properties is more complicated and logically involved than the ones for ( 1 ) and the KMRE. Characteristics of the predator equation PDE (5b) admits characteristics of the form d dt u ( t , τ 0 + t ) = - μ ( x , τ 0 + t ) u ( t , τ 0 + t ) , t ≤ τ , τ : = τ 0 + t , d d τ u ( τ + t 0 , τ ) = - μ ( x , τ ) u ( τ + t 0 , τ ) , t > τ , t : = τ + t 0 , 6 where τ 0 is the initial age and t 0 is the birth time (see Fig. 1 ). Fig. 1. Open in a new tab The age-time diagram for the model (7a)–(7h). The initial age distribution of the predator population is supported on the interval [ 0 , τ 0 , max ] marked with a thick green line. The time T is the maximal time until which we intend to compute the solution. Therefore, the interval [ 0 , τ T , max ] , where τ T , max : = τ 0 , max + T , is the age interval where the predator population at time T is supported With this in mind, we rewrite the system (5a)–(5g) as an infinite family of ODEs with renewal: Discretization The discretization of the model (7b)–(7h) serves two purposes: the numerical implementation and the proof of the existence of the solution. The age line is discretized with step h . The prey ODE (5a) and the family of predator ODEs ( 6 ) are discretized in time using the Forward Euler method with the same step h . This guarantees the propagation of generations of predators along the corresponding characteristics. The integral (5c) for the number of newborn predators and the integrals for the numbers of the juvenile and adult predators (5d) are calculated using the trapezoidal rule. We aim to compute a numerical solution till the maximal time T . Since the initial age distribution has a compact support [ 0 , τ 0 , max ] , the maximal predator age at time T is τ T , max : = τ 0 , max + T . We assume that the maturation age τ ∗ , the maximal initial age τ 0 , max , and the maximal time T are h -commensurable, i.e., τ ∗ = N 1 h , τ 0 , max = N 0 h , and T = N h for positive integers N , N 0 , and N 1 . Then the maximal age at time T is τ T , max = N 2 h where N 2 = N 0 + N . The numerical solution to the prey equation at time nh is denoted by X [ n ]. The numerical solution to the predator equation at time nh and age kh is denoted by U [ n , k ]. The resulting numerical scheme is: X [ n + 1 ] = X [ n ] 1 + h r - a X [ n ] + s Y 1 [ n ] - b Y 2 [ n ] , 8 U [ n + 1 , k ] = U [ n , k - 1 ] 1 - h μ ( X [ n ] , h ( k - 1 ) , k = 1 , … , N 2 , 9 U [ n + 1 , 0 ] = h ( 1 2 B ( X [ n ] , 0 ) U [ n , 0 ] + 1 2 B ( X [ n ] , τ T , max ) U [ n , N 2 ] + ∑ k = 1 N 2 - 1 B ( X [ n ] , h k ) U [ n , k ] ) ) , 10 Y 1 [ n + 1 ] = h 1 2 U [ n , 0 ] + 1 2 U [ n , N 1 ] + ∑ k = 1 N 1 - 1 U [ n , k ] ) , 11 Y 2 [ n + 1 ] = h 1 2 U [ n , N 1 ] + 1 2 U [ n , N 2 ] + ∑ k = N 1 + 1 N 2 - 1 U [ n , k ] ) . 12 The initial conditions for this scheme are X [ 0 ] = x 0 , 13 U [ 0 , k ] = u 0 ( h k ) , k = 0 , … , N 2 , 14 Y 1 [ 0 ] = h 1 2 u 0 ( 0 ) + 1 2 u 0 ( h N 1 ) + ∑ k = 1 N 1 - 1 u 0 ( h k ) , 15 Y 2 [ 0 ] = h 1 2 u 0 ( h N 1 ) + 1 2 u 0 ( h N 2 ) + ∑ k = N 1 + 1 N 2 - 1 u 0 ( h k ) . 16 The main theorem The goal of this section is to establish the existence, uniqueness, and positivity of the solution to the initial-value problem (7a)–(7h). We define a space X for pairs ( x , u ), X = { ( x , u ) | x ∈ R , u ∈ L 1 ( [ 0 , τ T max ] ) } , 17 with the norm ‖ ( x , u ) ‖ = | x | + ∫ 0 T , max | u | d τ . 18 Hence, X with the norm ( A3 ) is Banach because it is a direct product of Banach spaces R and L 1 ( [ 0 , τ T , max ] ) . Theorem 1 Let the initial prey population size x 0 be positive. Let the initial predator age density u 0 ( x ) be nonnegative, piecewise smooth, and compactly supported on [ 0 , τ 0 , max ] . Let the constants r , a , s , b , and μ M be positive and the constants k and g be nonnegative. Furthermore, let the function B ~ ( τ ) be continuous, positive, and bounded, and the function μ B ( τ ) be continuous, positive, and bounded at any finite τ . Let 0 ≤ t ≤ T . Then there exists a time interval [ 0 , T ∗ ] , T ∗ ≤ T , where the solution to (7a)–(7h) exists and is unique and positive. Moreover, the numerical solution by the scheme ( 8 )–( 16 ) converges to the exact solution to system (7a)–(7h) as the time and age step h → 0 , and the numerical error of the age distribution decays as O ( h ) at any fixed time t in the sense of the norm defined in ( 18 ). The proof of Theorem 1 is long and tedious. We sketch it now and elaborate it in Appendix A . The proof of the existence of the solution relies on the convergence of the numerical solution to a limit solution as the step h tends to zero. We show that the right-hand side of (7a)–(7d) is Lipschitz with respect to the norm ( A3 ) in a bounded region Ω R = { ( x , u ) ∈ X | ‖ ( x , u ) ‖ ≤ R } , 19 where R is a positive constant. Next, we fix step h and compute the numerical solution ( X h , U h ) till either time T ∗ , which is the minimum of T , and the exit time from the region Ω R . In addition, we compute a sequence of numerical solutions ( X h 2 - p , U h 2 - p ) till time t = n h ≤ T ∗ with steps h 2 - p , p = 1 , 2 , … and prove that ‖ ( X h 2 - p [ 2 p n ] , U h 2 - p [ 2 p n , · ] ) - ( X h [ n ] , U h [ n , · ] ) ‖ = O ( h ) , p = 1 , 2 , … . 20 Hence, this sequence of numerical solutions at time t = n h is Cauchy. Therefore, it converges to an element ( x ( t ) , u ( t ) ) ∈ X . The uniqueness follows from the Lipschitz continuity of the right-hand side of (7a)–(7d) and Grönwall’s inequality. The proof of positivity is conducted in four stages. First, we prove that x ( t ) is positive on [ 0 , T ∗ ] . Second, we argue that u ( t , τ ) is positive for t < τ . Third, we show that u ( t , 0) is positive. Finally, we conclude that u ( t , τ ) for t > τ is positive. Remark 1 The x -component and u along the characteristics are continuously differentiable as immediately follows from time stepping. The function u ( t , τ ) has a discontinuity along the characteristic t = τ unless u 0 ( 0 ) = ∫ 0 ∞ B ( x ( 0 ) , τ ) u 0 ( τ ) d τ . If u 0 ( τ ) has discontinuities at τ 1 , … , τ S , u ( t , τ ) also has discontinuities propagating along the characteristics t = τ - τ j , 1 ≤ j ≤ S . Therefore, the solution is well-defined as the solution to (7a)–(7h) everywhere except for the lines t = τ and t = τ - τ j , 1 ≤ j ≤ S . Reduction to an ODE model The goal of this section is to reduce the model (5a)–(5g) to an ODE model and compare the resulting ODE with ( 1 ) (Li et al. 2022 ). ODEs governing the population sizes of juvenile and adult predators are obtained by integrating PDE (5b) in age τ over the intervals [ 0 , τ ∗ ) and [ τ ∗ , ∞ ) respectively: d dt ∫ 0 ∗ u ( t , τ ) d τ + u ( t , τ ∗ ) - u ( t , 0 ) = - ∫ 0 ∗ μ ( x , τ ) u ( t , τ ) d τ , 21 d dt ∫ τ ∗ ∞ u ( t , τ ) d τ - u ( t , τ ∗ ) = - ∫ τ ∗ ∞ μ ( x , τ ) u ( t , τ ) d τ . 22 The integrals in the right-hand sides of ( 21 ) and ( 22 ) are expanded using the definition (5f) of the death rate function. In the calculations below, we replace φ { τ < τ ∗ } ( τ ) with the indicator function of [ 0 , τ ∗ ) for the sake of simplicity: ∫ 0 ∗ μ ( x , τ ) u ( t , τ ) d τ = g x y 1 + ∫ 0 ∗ μ B ( τ ) u ( t , τ ) d τ + μ M e - ρ x y 1 , 23 ∫ τ ∗ ∞ μ ( x , τ ) u ( t , τ ) d τ = ∫ τ ∗ ∞ μ B ( τ ) u ( t , τ ) d τ + μ M e - ρ x y 2 . 24 The expression for u ( t , 0) is obtained using the definition of the birth rate function (5e) and the assumption that B ~ ( τ ) is zero on [ 0 , τ ∗ ) : u ( t , 0 ) = ∫ 0 ∞ B ( x , τ ) u ( t , τ ) d τ = k x y 2 + ( 1 - e - ζ x ) ∫ τ ∗ ∞ B ~ ( τ ) u ( t , τ ) d τ . 25 Recalling the definitions of y 1 and y 2 , (5d), and using ( 21 )–( 25 ) we obtain the following ODEs for the juvenile, y 1 , and adult, y 2 , predator population sizes y 1 ′ + u ( t , τ ∗ ) - k x y 2 - ( 1 - e - ζ x ) ∫ τ ∗ ∞ B ~ ( τ ) u ( t , τ ) d τ = - g x y 1 - ∫ 0 ∗ μ B ( τ ) u ( t , τ ) d τ - μ M e - ρ x y 1 , 26 y 2 ′ - u ( t , τ ∗ ) = - ∫ τ ∗ ∞ μ B ( τ ) u ( t , τ ) d τ - μ M e - ρ x y 2 . 27 The term u ( t , τ ∗ ) is the number of juvenile predators becoming adults per time unit. Therefore, we define the transition rate from juvenile to adult predator as D = u ¯ ( τ ∗ ) y 1 , 28 where u ¯ ( τ ∗ ) is the value of u ( t , τ ∗ ) at the equilibrium. To eliminate integrals in ( 26 ) and ( 27 ), we introduce age-averaged birth rate b 2 for adult predators and age-averaged death rates m 1 and m 2 for juvenile and adult predators respectively: b 2 : = ∫ τ ∗ ∞ B ~ ( τ ) u ( t , τ ) d τ ∫ τ ∗ ∞ u ( t , τ ) d τ , 29 m 1 : = ∫ 0 ∗ μ B ( τ ) u ( t , τ ) d τ ∫ 0 ∗ u ( t , τ ) d τ , 30 m 2 : = ∫ τ ∗ ∞ μ B ( τ ) u ( t , τ ) d τ ∫ τ ∗ ∞ u ( t , τ ) d τ . 31 The resulting ODE system becomes ODE model (32a)–(32c) exposes to connection with ODE ( 1 )  Li et al. ( 2022 ). To obtain ( 1 ) from (32a)–(32c), one needs to remove the natural birth rate term proportional to b 2 y 2 and the hunger death rate terms proportional to μ M e - ρ x . In turn, to obtain the ODE system (32a)–(32c) from the proposed model (5a)–(5g), we replaced the birth and death rate functions with their age averages and replaced the smoothed indicator functions with the sharp ones. Reduction to a DDE model We also can reduce the model (5a)–(5g) to a delayed differential equation (DDE). We proceed as in Section 4 except for a different approximation of the term u ( t , τ ∗ ) in ( 26 )–( 27 ). Since the death rate μ ( x , τ ) depends on x and x depends on u , we cannot integrate the predator PDE (5b) along the characteristics τ = t + τ 0 and τ = t - t 0 explicitly. Note that the death rates in Mishra et al. ( 2024 ) and Mohr et al. ( 2014 ) were assumed independent of the prey population size – that’s why the integration was readily done in these works. Instead, we will apply the trapezoidal rule in x to obtain u ( t , τ ∗ ) by integrating along characteristics. Let t ≥ τ ∗ . Then the characteristic along which we will integrate can be parametrized as t - τ ∗ + s , 0 ≤ s ≤ τ ∗ , where t - τ ∗ is the birth time of the corresponding generation. Dividing (7c) by u ( t , τ ) and replacing the smooth indicator function with the sharp one we get: d ds log u ( t - τ ∗ + s , s ) = - g x ( t - τ ∗ + s ) - μ B ( s ) - μ M e - ρ x ( t - τ ∗ + s ) . 33 Integrating ( 33 ) from 0 to τ ∗ in s using the trapezoidal rule for x we obtain: u ( t , τ ∗ ) = u ( t - τ ∗ , 0 ) exp - g τ ∗ 2 [ x ( t - τ ∗ ) + x ( t ) ] - M B - μ M τ ∗ 2 e - ρ x ( t - τ ∗ ) + e - ρ x ( t ) . 34 The initial density u ( t - τ ∗ , 0 ) is found as in ( 25 ) except for t is replaced with t - τ ∗ . The symbol M B denotes the total death rate for juvenile predators: M B : = ∫ 0 τ μ B ( τ ) d τ . 35 We obtain u ( t , τ ∗ ) in the case t < τ ∗ in a similar manner. Thus, the DDE model is given by where the constants b 2 , m 1 and m 2 are defined by ( 29 ), ( 30 ), and ( 31 ) respectively. In contrast to ODE model (32a)–(32c), the DDE model (36a)–(36e) explicitly contains the maturation age parameter τ ∗ . Furthermore, this parameter enters the system of equations for the equilibrium of (36a)–(36e), where we set the right-hand sides to zero and x , y 1 , and y 2 to their equilibrium values: note the argument of the exponent in the expression for u ( t , τ ∗ ) for t ≥ τ ∗ . Therefore, the equilibria of the ODE and DDE models will not be exactly the same. Numerics We discretize the proposed model in the form (7a)–(7h) as described in Section 3.2 . The resulting system is ( 8 )–( 12 ) with initial conditions ( 14 )–( 16 ). The integration of the ODE model (32a)–(32c) and the DDE model (36a)–(36e) have been done using Matlab’s ode45 and dde23 respectively. Our codes are available on GitHub (Suarez and Cameron 2025 ). Settings For the sake of computational convenience, we introduce the maximal predator lifespan L = 30 . This is equivalent to having the base death rate function B ~ ( τ ) equal to infinity for τ > L . In Section 2 , we did not specify the base birth rate function B ~ ( τ ) and the base death rate function μ B ( τ ) to keep the model more general. Now, we define these functions as B ~ ( τ ) = 0 , τ < τ ∗ b p ( e - b ep ( τ - τ ∗ ) + 1 ) , τ ≥ τ ∗ , μ B ( τ ) = d p e d ep ( τ - L ) . 37 This lifespan L = 30 is large enough so that the predator’s age density is negligibly small at τ = L . The initial condition for the prey is set to x ( 0 ) = 0.5 . 38 The initial age density for the predator is chosen to be u 0 ( τ ) = 0.1 , τ ∈ [ 0 , τ ∗ ) 0.05 , τ ∈ [ τ ∗ , L ] 0 , otherwise . 39 These settings leave us with 15 parameters whose biological meaning and values are listed in Table 1 . The selected parameter values, when applicable, are borrowed from  Li et al. ( 2022 ) and are not attached to any particular predator-prey system. The value of b , the consumption rate of prey by adult predators, is increased from 0.4 to 0.8. The ranges for the parameters are chosen around those selected values to keep the order of magnitude of the selected value. The exception is the smoothness parameter of the indicator function ν . When ν = 100 , the indicator function is almost the Heaviside function. When ν = 1 , it changes gradually. For example, φ { τ ≥ 1 } ( 0 ) takes values of approximately 0.27 and 3.7 · 10 - 44 for ν = 1 and 100, respectively. Finding stable and unstable equilibria Stable and unstable equilibria of ( 8 )–( 12 ) are found using Newton’s method. We define a vector-function F ( X ∗ , U ∗ ) , F : R N 2 + 2 → R N 2 + 2 by subtracting the right-hand sides of ( 8 ) and ( 9 ) from their left-hand sides and replacing X [ · ] with X ∗ and U [ · , k ] with U ∗ [ k ] . Then, we set F ( X ∗ , U ∗ ) = 0 and solve this equation by Newton’s iteration. The ( N 2 + 2 ) × ( N 2 + 2 ) Jacobian matrix J is approximated by forward differences J : , 0 = 1 ϵ F ( X ∗ + ϵ , U ∗ ) - F ( X ∗ , U ∗ ) , 40 J : , k + 1 = 1 ϵ F ( X ∗ , U ∗ + e k ϵ U ∗ [ k ] ) - F ( X ∗ , U ∗ ) , k = 0 , … , N 2 41 where ϵ = 10 - 6 , “:" denotes all entries from 0 to N 2 + 1 , and e k is the standard k th basis vector in R N 2 + 1 . The warm start is obtained by a long integration of ( 8 )–( 12 ) at a set of parameter values where the equilibrium is stable and then by a continuation along a parameter under investigation. The stability of the equilibrium was assessed by computing the eigenvalues of the Jacobian matrix J and checking that all of them had negative real parts. Equilibria of the ODE and DDE systems were found using Matlab’s lsqnonlin, a nonlinear least-squares solver. The stability of the equilibria of the ODE model was checked by computing the eigenvalues of the Jacobian. The stability of the equilibria of the DDE model was checked by perturbing the equilibrium and examining where the system settled after a long time of integration. Finding limit cycles Periodic coexistence attractors exist at parameter values where the coexistence equilibrium is unstable. Periodic attractors, or limit cycles, are fixed points of Poincaré maps on appropriately chosen Poincaré sections. As a Poincaré section for the discretized age-structured model ( 8 )–( 12 ), we choose the hyperplane in R N 2 + 2 X = X ∗ , where X ∗ is the first component of the unstable coexistence equilibrium. The Poincaré map G : R N 2 + 1 → R N 2 + 1 maps the point where a trajectory of ( 8 )–( 12 ) crosses the Poincaré section to the next crossing in the same direction. The point of intersection is found using time integration of ( 8 )–( 12 ) and linear interpolation. The fixed points of the Poincaré map are found using the Levenberg-Marquardt nonlinear solver (Nocedal and Wright 2006 , Section 10.3). The objective function for the Levenberg-Marquardt method is defined as f ( U ) = 1 2 ‖ U - G ( U ) ‖ 2 2 . The Jacobian of U - G ( U ) is found using forward differences similar to ( 41 ). The Levenberg-Marquardt method reliably found limit cycles, while Newton’s method failed to do so in our settings because it required warmer starts than it was practical for us to provide. The limit cycles for the ODE system were found similarly. Long-time integration found the limit cycles for the DDE system. Results We subjected the age-structured model (5b)–(5g) to a detailed numerical investigation. First, we performed the Latin Hypercube Sampling (Audze and Eglãjs 1977 ; McKay et al. 1979 ; Iman et al. 1981 ) in the 15-dimensional parameter space. Then we conducted the Linear Discriminant Analysis (LDA) (Duda et al. 2001 ) and found which parameters affect the type of attractor the most. Next, we carried out a thorough investigation of the dynamical behavior of the age-structured model (5b)–(5g) depending on the maturation age τ ∗ and the consumption rate of the juvenile predator by the prey g , the two most important parameters. Finally, we repeated the investigation of the dynamical behavior of the corresponding ODE, (32a)–(32c), and DDE, (36a)–(36e), models in the same region of the ( τ ∗ , g ) -plane and compared the calculated phase diagrams for all three models. Latin Hypercube Sampling The Latin Hypercube Sampling aims to examine the types of dynamical behavior that may occur in the system depending on its parameters. As we have shown in Section 3.3 , the solution to (5b)–(5g) exists, is unique, and is positive. Latin Hypercube Samples were generated in 15-dimensional space formed by the direct product of the intervals specified in Table 1 . The initial condition for each parameter set was the same as in ( 38 ) and ( 39 ). At each parameter set, the numerical integration continued either till time T max = 500 or till any of the solution components exceeded 1000. In the latter case, we declared a blow-up. The choice of the time step for Latin Hypercube is a trade-off between reliability and runtime. We monitored the positivity of the solution components throughout the run. Negative values of u or x result from numerical errors and mean that one needs to reduce the time step. By such trial and error, we chose the time step Δ t = 0.005 . We conducted two runs of Latin Hypercube Sampling of 10,000 samples each. Four types of long-term behavior were observed: These types of behavior are illustrated and a summarising bar graph is presented in Fig. 2 . Fig. 2. Open in a new tab The results of Latin Hypercube Sampling. The intervals for the parameters are specified in Table 1 The Coexistence Equilibrium, Coexistence Periodic, and Predator-Free Attractors are biologically relevant types of behavior, while the blow-up is not. The blow-up scenarios are characterized by oscillations of the prey and predator population sizes of increasing amplitude. Such growing oscillations seem to be possible due to the unbounded predator birth-rate function B ( x , t ) in equation (5e), the unbounded prey birth rate r + s y 1 in the prey ODE (5a), and time delay effects. A positive feedback loop may arise in the model when the consumption rate of the juvenile predator, g , is small, the maturation age τ ∗ is large, and the growth rate of the prey due to eating juvenile predators, s , is large. Then, the more prey there are, the more prey is eaten by the adult predators, the more the juvenile predators are born, and the more prey grows due to feeding on juvenile predators without damaging them much. More biologically plausible birth rates must contain saturation, as no species can reproduce infinitely fast. For example, one may consider the following modifications to the prey ODE and the predator birth rate function, (5a) and (5e), respectively: x ′ = x r - a x + s y ^ 1 tanh y 1 y ^ 1 - b y 2 , 42 B ( x , τ ) = k x ^ tanh x x ^ φ τ ≥ τ ∗ ( τ ) + B ~ ( 1 - e - ζ x ) , 43 where y ^ 1 and x ^ are parameters. The function tanh , a popular activation function in neural network-based smooth solution models to PDEs used, e.g., in Li et al. ( 2019 ), vanishes at zero, is approximately equal to its argument on [0, 1], and nearly reaches its upper bound when its argument is greater than 4. The resulting prey birth rate term, , is approximately equal to s y 1 when y 1 ≤ y ^ 1 , and is bounded from above by s y ^ 1 , when y 1 → ∞ . The predator birth rate term behaves likewise. We set x ^ = 20 and y ^ 1 = 10 and repeated Latin Hypercube Sampling in the 15D parameter space with parameter ranges from Table 1 . The system with saturated birth rates settled at Equilibrial Coexistence Attractor, Periodic Coexistence Attractor, and Predator-Free Attractor in 1631, 2529, and 5840 cases, respectively. No blow-up was registered. The existence, uniqueness and positivity of the age-structured model with modifications ( 42 ) and ( 43 ) can be proven as it is done for our original age-structured model (7a)–(7h). The proof of the boundedness of the solution to this modified age-structured model can be outlined as follows. The boundedness of the prey population follows from the comparison principle for ODEs and the dominance of the right-hand side of ( 42 ) by the function x ( r + s y 1 ^ ) - a x 2 for all positive x . Hence, the solution x ( t ) is bounded by a - 1 ( r + s y ^ 1 ) . The predator birth rate function ( 43 ) is bounded by B ¯ : = k x ^ + B from above. The predator death rate is bounded from below by μ ¯ : = d p exp ( - d ep L ) in our model. Therefore, the predator population will be dominated by the solution to the KMRE and hence will be bounded. We leave further improvements of the age-structured model and a rigorous proof of its properties for future work. Linear Discriminant Analysis Linear Discriminant Analysis (LDA) (or Multiple Discriminant Analysis (MDA)) is a classical linear supervised learning tool (Duda et al. 2001 ). It aims at finding a low-dimensional subspace such that data from different categories projected onto this subspace are separated the most while the projected data from the same categories are clustered the most. A description of LDA is found in Appendix B . We use LDA to discriminate parameter sets with four different types of long-term behavior and find which parameters affect the type of long-term behavior the most. Since there are four categories, the maximal dimension of the optimal subspace is three. We project the parameters onto an optimal plane spanned by the two dominant eigenvectors of the generalized eigenvalue problem described in Appendix B – see Fig. 3 . Fig. 3. Open in a new tab LDA projection from the 15D parameter space onto a 2D optimal subspace where the parameter sets leading to four different types of long-term behavior (Blow-Up, Coexistence Equilibrium, Coexistence periodic, and Predator-free attractor) are separated the most Note that the eigenvectors of the LDA generalized eigenvalue problem are not orthogonal. We orthonormalized them by one step of the Gram-Schmidt procedure. To understand which parameters affect the type of long-term behavior the most, we do the following calculation. Let X be n × d matrix whose rows are the parameter sets generated by the Latin Hypercube Sampling, n = 10 , 000 , and d = 15 in our case. Let R be a diagonal matrix whose diagonal entries R i are the differences between the maximal and minimal values of parameter i , i = 1 , … , d . Then, the matrix X can be decomposed as X = Z R + 1 n × 1 x min , 44 where x min is a row vector whose entries are the lower bounds for the parameter ranges, and Z is an n × d matrix whose entries take values between zero and 1. Therefore, to account for ranges, we take the d × 2 LDA projection matrix W with orthonormalized columns and multiply it by R on the left. Then, we plot bar graphs of RW and display the result in Fig. 4 . Fig. 4. Open in a new tab LDA range-adjusted bar graphs quantifying the importance of the model parameters for the type of long-term behavior This graph suggests that the first LDA is the most influenced, in decreasing order, by g , the consumption rate of juvenile predators by the prey, τ ∗ , predator maturation age, a , intraspecific competition rate of the prey, and r , the prey birth rate. The second LDA is the most affected, in decreasing order, by g , k , the reproduction rate of the predator, s , the consumption rate of prey by adult predators, and b , the rate of predation by the predator. Phase diagrams in the ( τ ∗ , g ) -plane Fig. 4 suggests that the two most influential parameters of our model on the long-term behavior type are the predator maturation age, τ ∗ , and the consumption rate of juvenile predators by the prey, g . We conducted a detailed investigation of the long-term behavior of the system (5a)–(5g) at ( τ ∗ , g ) ∈ [ 0 , 2 ] × [ 0 , 1 ] , and the rest of the parameters fixed at the selected values specified in the third column of Table 1 . We used a time/age step of h = 0.0125 . For each value of τ ∗ from 0.1 to 2 with step 0.05, we moved along the parameter g with step 0.01 from g = 1 down to g = 0 . At each ( τ ∗ , g ) , we first ran a simulation on the time interval 0 ≤ t ≤ 500 starting from the initial condition x ( 0 ) = 0.5 , u ( 0 , τ ) = 0.1 , τ ≤ τ ∗ , 0.05 , τ > τ ∗ . If this simulation resulted in the extinction of the predator, we declared the Predator-Free Attractor. Otherwise, we used the prey population size and the predator age density at time t = 500 as the initial guess for finding the equilibrium. We found the equilibrium as described in Section 6.2 . If the equilibrium was unstable, i.e., if the Jacobian ( 40 )–( 41 ) had an eigenvalue with a positive real part, we found the Periodic Coexistence attractor as described in Section 6.3 . We did so for two values of the indicator function smoothness parameter ν : ν = 100 corresponding to a sharp transition from juvenile to adult, and ν = 1 , describing a gradual transition. The resulting phase diagrams in the plane ( τ ∗ , g ) , bifurcation diagrams at τ ∗ = 1 , and ensembles of bifurcation diagrams in g corresponding to the grid values of τ ∗ are displayed in Fig. 5 . Comparing these diagrams at ν = 100 and ν = 1 , we observe that the gradual transition from juvenile to adult somewhat increases the region of the Predator-Free Attractor, slightly reshapes the region of the Periodic Coexistence Attractor, and slightly increases the amplitude of stable oscillations at large τ ∗ compared with those in the case of the sharp transition. We also observe an oscillatory feature in the upper branches of the bifurcation diagrams for the juvenile predator at both ν = 100 and ν = 1 . These features persisted throughout our time- and age-step refinement. We leave an investigation into this phenomenon for future work. Fig. 5. Open in a new tab Phase and bifurcation diagrams for the proposed model (5a)–(5g). A sharp (left) and a gradual (right) transition from juvenile to adult predator: ν = 100 and ν = 1 in (5g), respectively. Top row: Phase diagrams. Middle Row: Bifurcation diagrams at the maturation age τ ∗ = 1 . Bottom row: Ensembles of the bifurcation diagrams Age density We examined predator age-density at several pairs of maturation age τ ∗ and the consumption rate of juvenile predators by prey g at which the system admits the Equilibrial Coexistence Attractor. The remaining parameters were set to their selected values in Table 1 . The results are shown in Fig. 6 . The predator age-density decays rapidly with age and approaches zero at the age cutoff L = 30 , as desired. Fig. 6. Open in a new tab The predator age density at several pairs of values ( τ ∗ , g ) leading to the Equilibrial Coexistence Attractor. The rest of the parameters are at their selected values from Table 1 , and ν = 100 We have also extracted the predator age density at ( τ ∗ = 1 , g = 0.1 ) where the system settles on the Periodic Coexistence Attractor – see Fig. 7 . Fig. 7. Open in a new tab The predator age density at ( τ ∗ = 1 , g = 0.1 ) leading to the Periodic Coexistence Attractor. The rest of the parameters are at their selected values from Table 1 , and ν = 100 We observe rapidly decaying age density waves and the following alternation of peaks and minima of the prey and juvenile and adult predator population sizes: max y 2 → min x → min y 1 → min y 2 → max x → max y 1 → max y 2 . Comparison with the ODE and DDE models The proposed model (5a)–(5g) was reduced to ODE and DDE models, (32a)–(32c) and (36a)–(36e), respectively, by replacing the smooth indicator functions with the sharp ones, integration over age, and age-averaging (see Section 4 ). The trapezoidal rule was used to approximate the transition term from juvenile to adult predator in the DDE model, yielding delayed terms by the maturation age τ ∗ . The goal of this section is to numerically investigate the ODE and DDE models and compare their long-term behavior to that of the age-structured model (5a)–(5g). We set all parameters of (5a)–(5g), except for τ ∗ and g , to their selected values from Table 1 . The parameters ( τ ∗ , g ) ran through all values from the rectangle [ 0 , 2 ] × [ 0 , 1 ] . The parameters for the ODE system (32a)–(32c), D , the rate from juvenile to adult predator, ( 28 ), b 2 , the predator birth rate, ( 29 ), m 1 , the death rate for juvenile predators, ( 30 ), and m 2 , the death rate for adult predators, ( 31 ), were computed by age-averaging at the equilibria, stable or unstable, admitted by the system (5a)–(5g) with ν = 100 at each pair ( τ ∗ , g ) . The age-averaged values of D , b 2 , m 1 , and m 2 as functions of τ ∗ and g are displayed in Fig. 8 . Their ranges are D : [ 0.353 , 9.91 ] , b 2 : [ 0.0794 , 0.0835 ] , m 1 : [ 0.0200 , 0.0219 ] , m 2 : [ 0.0363 , 0.0546 ] . The DDE model (36a)–(36e) uses the same parameter values b 2 , m 1 , and m 2 as the ODE model. It does not involve the parameter D . Fig. 8. Open in a new tab The dependence of the age-averaged parameters of the ODE model (32a)–(32c) on the maturation age τ ∗ , and the consumption rate of juvenile predators by the prey, g . The white curves are the boundaries between the regions with the Periodic and Equilibrial Coexistence Attractors of the age-structured model Thus, with the selected parameter values from Table 1 and Fig. 8 , the ODE and DDE systems mimic the age-structured system (5a)–(5g) as closely as possible. The phase diagrams in ( τ ∗ , g ) for the ODE system (32a)–(32c) and the DDE system (36a)–(36e) are shown in Fig. 9 . Fig. 9. Open in a new tab (a): The phase diagram in the ( τ ∗ , g ) plane and a bifurcation diagram at τ ∗ = 2 for the ODE model (32a)–(32c) with parameters D , b 2 , m 1 , and m 2 age-averaged over the corresponding attractors of the proposed model (5a)–(5g) and displayed in Fig. 8 . (b): The phase diagram in the ( τ ∗ , g ) plane and a bifurcation diagram at τ ∗ = 1.3 for the DDE model (36a)–(36e) where parameters b 2 , m 1 , and m 2 age-averaged over the corresponding attractors of the proposed model (5a)–(5g) and displayed in Fig. 8 As in the age-structured model, the phase diagrams for the ODE and DDE systems have three regions corresponding to the Equilibrial Coexistence Attractor, the Periodic Coexistence Attractor, and the Predator-Free Attractor. The regions of the Predator-Free Attractor for the ODE and DDE models did not increase compared to the age-structured model (5a)–(5g) in Fig. 5 (top left). We remark that we did not run the ODE and DDE models in the Predator-Free region of the age-structured model with ν = 100 because we did not compute the parameters D , b 2 , m 1 , and m 2 in it. The region of the Periodic Coexistence Attractor is considerably smaller for the ODE system. This region for the DDE system is intermediate in size between those of the ODE and age-structured models. Fig. 10 superimposes the phase diagrams of the age-structured model (5a)–(5g) with the indicator function smoothness parameters ν = 1 and 100, the ODE system (32a)–(32c), and the DDE system (36a)–(36e). The age-structured model is much more prone to developing oscillations than the closely mimicking ODE and DDE systems. Fig. 10. Open in a new tab A comparison of phase diagrams on the ( τ ∗ , g ) plane of the age-structured model (5a)–(5g) with the indicator function smoothness coefficients ν = 100 and 1, the ODE model (32a)–(32c), and the DDE model (36a)–(36e) with parameters mimicking the age-structured model with ν = 100 Discussion We developed a predator-prey model (5a)–(5g) with an age-structured predator population and role reversal, proved the existence, uniqueness, and positivity of the initial-value problem for it, and investigated its long-term behavior numerically. In addition, we derived ODE and DDE models from it and compared their long-term behavior to that of the age-structured model with the corresponding parameter values. This study taught us several important lessons. The proposed model involves 15 parameters. Latin Hypercube Sampling has shown that, depending on these parameters, the system settles on the Predator-Free Attractor, Equilibrial Coexistence Attractor, and Periodic Coexistence Attractor in 22%, 19%, and 55% of cases, respectively, or blows up in approximately 4% of cases. The blow-up is not biologically relevant. It is enabled due to the terms k x φ { τ ≥ τ ∗ } ( τ ) in the birth rate function of the predator, (5e) and s x y 1 in the right-hand side of the prey ODE (5a). The replacement of the birth-rate function and the prey ODE with their saturated versions ( 42 ) and ( 43 ) resulted in the absence of blow-ups in the Latin Hypercube Sampling. Nonetheless, we left this caveat in our model for educational purposes. Our future models will necessarily involve saturation, preventing blow-ups. We performed a Linear Discriminant Analysis to find which parameters affect the type of long-term behavior the most. These are the maturation age τ ∗ and the consumption rate of the juvenile predator by the prey g . The importance of these parameters agrees with the main claim of Werner and Gilliam ( 1984 ): the maturation age of the predator and its ontogenetic niche shift in diet and interspecific interactions are significant determinants of the dynamics of simplified ecological systems. We fixed all parameters but τ ∗ and g at their selected values (see Table 1 ) and further investigated the system’s dynamics. Blow-up does not occur in these settings, and the system approaches one of the three attractors depending on τ ∗ and g . The phase diagrams in Fig. 5 in the ( τ ∗ , g ) -plane reveals the following. If the maturation age is small enough, the system reaches a coexistence equilibrium at all values of g ∈ [ 0 , 1 ] . As the maturation age τ ∗ increases, the predator becomes extinct at large g , and stable oscillations develop at small g . The region of the Predator-Free Attractor occupies the top right corner of the phase diagrams where both τ ∗ and g are large, while the Periodic Coexistence Attractor occupies the bottom right corner where τ ∗ is large and g is small. This suggests that the long maturation of the predator and significant consumption of its juveniles create a juvenile bottleneck that is hard to overcome and results in the predator’s extinction. This result is in agreement with Werner and Gilliam ( 1984 ). The effect of a gradual rather than a sharp transition from juvenile to adult is a minor increase of the Predator-Free Attractor region, a minor reshaping of the Periodic Coexistence Attractor region, and a minor increase of the amplitude stable oscillations at large τ ∗ . Overall, this effect is mild. We extracted the predator age density at a collection of pairs ( τ ∗ , g ) leading to an Equilibrial Coexistence Attractor and at a pair ( τ ∗ , g ) admitting the Periodic Coexistence Attractor. The age densities at Equilibrial Coexistence Attractors monotonically decay, while they feature decaying waves at the periodic Coexistence Attractor. This behavior of the age density is biologically relevant. To gauge the effect of the age-structured predator population, we derived ODE and DDE models. Furthermore, we computed the parameter values for these models by age-averaging the solution to the age-structured model at its coexistence equilibrium, whether stable or unstable. The phase diagram for the resulting ODE model has a much smaller area of the Periodic Coexistence Attractor and a much larger region of the Equilibrial Coexistence Attractor. The sizes of these regions in the phase diagram of the DDE model are intermediate between those of the ODE and age-structured models. Therefore, the age structure of the predator population promotes oscillatory behavior. Numerical simulation of the ODE and DDE models is much faster than that of the age-structured model. Time and age step refinement in the age-structured model by a factor of two increases runtime by a factor of four. The computation time for the phase diagrams in Fig. 10, which involve finding equilibria and limit cycles, was several days for the age-structured model, compared with several hours for the ODE and DDE models. The finding that different models of the same system can yield qualitatively different predictions is not new. For example, Cantrell et al. ( 2012 ) showed that spatially explicit and spatially implicit models of a one-dimensional two-patch system yield different predictions for the survival of a population migrating between the two patches. Conclusion The proposed predator-prey model (5a)–(5g) with a role reversal and age-structured predator population captures the role of predators’ gradual maturation and the ontogenetic niche shifts in their diet and their interaction with the prey. The long-term behavior of this system at the default parameter values agrees with our knowledge about some real predator-prey systems. A comparison of the long-term behavior of the age-structured model with the mimicking ODE- and DDE-based models demonstrates that the model type strongly affects the type of coexistence attractor the system admits. Because the age-structured model includes elements of biological reality that are absent in the ODEs and DDEs, one should use caution when drawing conclusions from the latter model types. The proposed model contains a number of imperfections. First, it admits blow-up for a small fraction of parameter cases, and we have suggested some strategies for revising the model that should avoid this. Second, the model does not account for Allee effects (Allee et al. 1949 ), demographic stochasticity, and other ecological processes important at small population sizes. Consequently, the lower bounds of oscillatory solutions to the coupled predator-prey system may be very low, which is ecologically implausible. Solving this issue would require major changes to the model formulation. Nevertheless, we have learned a lot from our study of this model and leave further improvements to future work. Acknowledgements The work of M.C. was partially supported by the AFOSR MURI grant FA9550-20-1-0397. W.F.F. acknowledges support from the U.S. National Science Foundation (DMS2451241). The work of D.L. was partially supported by the Simons Foundation. Proof of Theorem 1 First of all, we remind the reader that the age-structured model (5a)–(5g) can be rewritten as a system of infinitely many ODEs, (7a)–(7h). The first ODE, (7a), describes the dynamics of the prey population size, while the two infinite families of ODEs, (7b) and (7c), describe the dynamics of generations of predators present at the initial time t = 0 and born at t > 0 , respectively. Compact support of the predator age density. Lemma 1 The age density of the predator is compactly supported at any finite time t . Proof The initial age density u 0 ( τ ) is bounded and has a compact support. The age density u ( t , τ ) propagates along characteristics (7b)–(7c). Equation (7b) implies that u ( t , t + τ 0 ) = 0 whenever u ( 0 , τ 0 ) = u 0 ( τ 0 ) = 0 . Hence, if u 0 ( τ ) = 0 on [ τ 0 , max , + ∞ ) then at a fixed t , u ( t , τ ) = 0 on [ τ 0 , max + t , + ∞ ) . □ For brevity, we denote the maximal possible age at time t by τ max ( t ) : τ max ( t ) : = τ 0 , max + t . A1 A Banach space for prey population size and predator age density. To analyze the solutions to the initial value problem (5a)–(5g), we define a space for the prey population size x and the predator age density u . Lemma 2 The space X = { ( x , u ) | x ∈ R , u ∈ L 1 ( [ 0 , τ max ( t ) ] ) } A2 with the norm ‖ ( x , u ) ‖ X = | x | + ∫ 0 τ max ( t ) | u | d τ . A3 is a Banach space at each fixed finite time t . Proof Indeed, L 1 ( [ 0 , τ max ( t ) ] ) is Banach and hence X is a direct product of Banach spaces – see (Adams and Fournier 2003 ). □ Lipschitz continuity. Lemma 3 The right-hand side F : = ( f 1 , f 2 ) of (7a)–(7c) is locally Lipschitz with respect to the norm ( A3 ), i. e. for any ( x , u ) , ( y , v ) such that ‖ ( x , u ) ‖ X ≤ R , ‖ ( y , v ) ‖ X ≤ R , A4 there exists a Lipschitz constant L F , R dependent on F and R but independent of ( x , u ) and ( y , v ) such that ‖ F ( x , u ) - F ( y , v ) ‖ X ≤ L F , R ‖ ( x , u ) - ( y , v ) ‖ X . A5 Proof We will compute | f 1 ( x , u ) - f 1 ( y , v ) | and ‖ f 2 ( x , u ) - f 2 ( y , v ) ‖ L 1 ( [ 0 , τ max ( t ) ] . A6 For brevity, we will use the notation ‖ f 2 ( x , u ) - f 2 ( y , v ) ‖ 1 : = ‖ f 2 ( x , u ) - f 2 ( y , v ) ‖ L 1 ( [ 0 , τ max ( t ) ] . A7 We start with f 1 , the right-hand side of (7a): | f 1 ( x , u ) - f 1 ( y , v ) | = = r ( x - y ) - a ( x 2 - y 2 ) + s x ∫ 0 ∗ u ( t , τ ) d τ - y ∫ 0 ∗ v ( t , τ ) d τ - b x ∫ τ ∗ τ max ( T ) u ( t , τ ) d τ - y ∫ τ ∗ τ max ( T ) v ( t , τ ) d τ ≤ | x - y | r + a ( x + y ) + s x ∫ 0 ∗ u ( t , τ ) d τ - y ∫ 0 ∗ v ( t , τ ) d τ + b x ∫ τ ∗ τ max ( T ) u ( t , τ ) d τ - y ∫ τ ∗ τ max ( T ) v ( t , τ ) d τ ≤ | x - y | r + a ( x + y ) + max { s , b } | x | ∫ 0 max ( t ) | u ( t , τ ) - v ( t , τ ) | d τ + max { s , b } | y - x | ∫ 0 τ max ( t ) | v ( t , τ ) | d τ ≤ | x - y | r + 2 a R + max { s , b } R | y - x | + ‖ u - v ‖ 1 . A8 We continue with f 2 ( x , u ) = - μ ( x , τ ) u , the right-hand side of (7b) and (7c), where the death rate function μ ( x , τ ) is given by (7f). The base death rate term in (7f), μ B ( τ ) , is continuous and bounded at every finite τ by the statement of Theorem 1 . Therefore, we assume that 0 ≤ μ B ( τ ) ≤ M B on 0 ≤ τ ≤ τ max ( t ) . Hence, ‖ f 2 ( x , u ) - f 2 ( y , v ) ‖ 1 = ∫ 0 τ max ( T ) | f 2 ( x , u ) - f 2 ( y , v ) | d τ ≤ g | x | ∫ 0 ∗ | u ( t , τ ) - v ( t , τ ) | d τ + g | x - y | ∫ 0 ∗ | u ( t , τ ) | d τ + M B ∫ 0 τ max ( t ) | u ( t , τ ) - v ( t , τ ) | d τ + μ M e - ρ x ∫ 0 τ max ( T ) | u ( t , τ ) - v ( t , τ ) | d τ + μ M e - ρ x - e - ρ y ∫ 0 τ max ( t ) | v ( t , τ ) | d τ ≤ g R ‖ u - v ‖ 1 + g R | x - y | + M B ‖ u - v ‖ 1 + μ M max { 1 , e - ρ R } ‖ u - v ‖ 1 + μ M R | x - y | W ρ , R , A9 where W ρ , R is the Lipschitz constant for e - ρ x in | x | ≤ R . Putting inequalities ( A8 ) and ( A9 ) together, we get ‖ F ( x , u ) - F ( y , v ) ‖ X = | f 1 ( x , u ) - f 1 ( y , v ) | + ‖ f 2 ( x , u ) - f 2 ( y , v ) ‖ 1 ≤ L F , R ‖ ( x , u ) - ( y , v ) ‖ X A10 where the Lipschitz constant L F , R is L F , R : = r + M B + R 2 a + max { s , b } + g + μ M W ρ , R + μ M max { 1 , e - ρ R } . A11 This completes the proof of Lemma 3 . □ How long does a numerical solution stay in a ball? The Lipschitz continuity of F , ( A5 ), and the fact that F ( 0 , 0 ) = ( 0 , 0 ) imply that ‖ F ( x , u ) ‖ X = | f 1 | + ‖ f 2 ‖ 1 ≤ L F , R ‖ ( x , u ) ‖ X ≤ L F , R R . A12 This allows us to guarantee that a numerical solution with any finite time step exists for at least a minimal time independent of the step size, if the step size is sufficiently small. Lemma 4 Let ( X , U ) be a numerical solution to ( 8 )–( 14 ) with a time step h . Let ‖ ( X , U ) ‖ X ≤ R 0 . Then the numerical solution will remain in the ball ‖ ( X , U ) ‖ X ≤ R where R > R 0 at least for some minimal positive time independent of h . Proof For the x -component of the solution we have: | X [ k + 1 ] | = | X [ k ] | + h | f 1 ( X [ k ] , U [ k ] ) | ≤ | X [ k ] | + h L F , R R ‖ ( X [ k ] , U [ k , · ] ) ‖ X . A13 Next, we estimate U [ k + 1 , 0 ] : U [ k + 1 , 0 ] = h Trap ( B ( X [ k ] , · ) U [ k , · ] ) ≤ B max ‖ U [ k , · ] ‖ 1 + O ( h ) ≤ ( B max + 1 ) ‖ ( X [ k ] , U [ k ] ) ‖ X , where h Trap ( · ) denotes the trapezoidal rule with step h applied to ( · ) and B max is the maximum of B ( X , τ ) in the ball ‖ ( X , U ) ‖ X ≤ R . The, we bound U [ k + 1 , j ] for j ≥ 1 : U [ k + 1 , j ] = U [ k , j - 1 ] + h f 2 ( X [ k ] , h ( j - 1 ) , U [ j - 1 ] ) . Integration over age, we get ‖ U [ k + 1 , · ] ‖ 1 ≤ h C 0 ‖ U [ k , · ] ‖ 1 + ‖ U [ k , · ] ‖ 1 + h L F , R R ‖ ( X [ k ] , U [ k ] ) ‖ X . A14 Adding ( A13 ) and ( A14 ) we obtain: ‖ ( X [ k + 1 ] , U [ k + 1 ] ) ‖ X ≤ ‖ ( X [ k ] , U [ k ] ) ‖ X ( 1 + C h ) , A15 where C : = 2 L F , R R + B max + 1 . Therefore, ‖ ( X [ k ] , U [ k ] ) ‖ X ≤ ‖ ( X [ 0 ] , U [ 0 ] ) ‖ X ( 1 + C h ) k ≤ R 0 e Ckh . A16 This means that the numerical solution ( X [ k ], U [ k ]) will remain in the ball ‖ ( X [ k ] , U [ k ] ) ‖ X ≤ R for the time at least T : = k h = 1 2 L F , R R + B max + 1 log R R 0 > 0 . A17 □ Existence. Lemma 5 Consider a time-space cylinder C : = { ( t , x , u ) | t ∈ [ 0 , T ] , ‖ ( x , u ) ‖ X ≤ R } . A18 Let the initial condition satisfy ‖ ( x 0 , u 0 ( τ ) ) ‖ X = R 0 < R . Then the initial value problem (7a)–(7h) has a solution on the time interval [ 0 , min { T , t R } ] where t R = inf { t ≥ 0 | ‖ ( x ( t ) , u ( t , τ ) ) ‖ X ≥ R } . A19 Proof We will construct a sequence of approximate, or numerical, solutions to (7a)–(7h) and show that this sequence is Cauchy and hence converges. This limit is the desired solution to the initial value problem (5a)–(5g). A numerical solution is obtained using the discretization ( 8 )–( 16 ) and time-stepping. We fix a large N , define h : = T / N , and construct a sequence of numerical solutions with time and age step sizes 2 - p h , where p = 0 , 1 , 2 , … . We denote the numerical solution with the time and age step 2 - p h by ( X p , U p ) . A continuous numerical solution in the age variable at any fixed time t = n h where n is small enough is defined by linear interpolation of ( X p [ n ] , U p [ n , · ] ) in age. We terminate time stepping as soon as the X -norm of the numerical solution with any time step 2 - p h , p = 0 , 1 , 2 , … , ( X 0 [ n ] , U 0 [ n , · ] ) , exceeds R , or as time reaches T = N h , whichever event happens first. Lemma 4 guarantees that the termination time is bounded from below. Let the terminal time be T 1 = N 1 h . Our goal is to prove that for any p ∈ N and any 1 ≤ n ≤ N 1 , ( X 0 [ n ] , U 0 [ n , · ] ) - ( X p [ 2 p n ] , U p [ 2 p n , · ] ) X ≤ A h , A20 where A is a constant independent of p . Note that 2 p n time steps of the numerical recurrence ( 8 )–( 16 ) with steps h and h 2 - p yield approximate solutions ( X 0 [ n ] , U 0 [ n , · ] ) and ( X p [ 2 p n ] , U p [ 2 p n , · ] ) to (7a)–(7h) at the same time t = n h . Step 1. The discrepancy between the coarse and fine mesh solutions over the first step. The first step toward this goal is to show that for any p ∈ N , ( X 0 [ 1 ] , U 0 [ 1 , · ] ) - ( X p [ 2 p ] , U p [ 2 p , · ] ) X ≤ A 1 h 2 , A21 where A 1 is a constant independent of p . The x -components of numerical solutions over the time interval h on the coarse h , and the fine, 2 - p h , meshes are, respectively, X 0 [ 1 ] = x 0 + h f 1 ( x 0 , u 0 ( · ) ) , X p [ 1 ] = x 0 + 2 - p h f 1 ( x 0 , u 0 ( · ) ) , ⋮ X p [ 2 p ] = X p [ 2 p - 1 ] + 2 - p h f 1 ( X p [ 2 p - 1 ] , U p [ 2 p - 1 , · ] ) . For the u -components of these solutions we have: U 0 [ 1 , k ] = u 0 ( h ( k - 1 ) ) + h f 2 ( x 0 , u 0 ( h ( k - 1 ) ) , h ( k - 1 ) ) , k ≥ 1 , U 0 [ 1 , 0 ] = h Trap ( B ( x 0 , · ) u 0 ( · ) ) , U p [ 1 , j ] = u 0 ( 2 p h ( j - 1 ) ) + 2 - p h f 2 ( x 0 , u 0 ( 2 p h ( j - 1 ) ) , 2 p h ( j - 1 ) ) , j ≥ 1 , U p [ 1 , 0 ] = 2 - p h Trap ( B ( x 0 , · ) u 0 ( · ) ) , ⋮ U p [ 2 p , j ] = U p [ 2 p - 1 , j - 1 ] + 2 - p h f 2 ( X p [ 2 p - 1 ] , U p [ 2 p - 1 , j - 1 ] , ( 2 p - 1 ) h ( j - 1 ) ) , j ≥ 1 , U p [ 2 p , 0 ] = 2 - p h Trap ( B ( X p [ 2 p - 1 ] , · ) U p [ 2 p - 1 , · ] ) . The discrepancy in x . Now we bound the difference between X 0 [ 1 ] and X p [ 2 p ] : | X p [ 2 p ] - X 0 [ 1 ] | ≤ 2 - p h ∑ l = 1 2 p - 1 | f 1 ( X p [ l ] , U p [ l , · ] ) - f 1 ( x 0 , u 0 ) | ≤ 2 - p h L F , R ∑ l = 1 2 p - 1 | X p [ l ] - x 0 | + ‖ U p [ l , · ] - u 0 ‖ 1 . A22 The Lipschitz continuity of F , ( A5 ), and the fact that F ( 0 , 0 ) = ( 0 , 0 ) imply that ‖ F ( x , u ) ‖ X = | f 1 | + ‖ f 2 ‖ 1 ≤ L F , R ‖ ( x , u ) ‖ X ≤ L F , R R . A23 The bound ( A23 ) allows us to bound | X p [ l ] - x 0 | , 0 ≤ l ≤ 2 p - 1 , as follows: | X p [ l ] - x 0 | = x 0 + 2 - p h ∑ i = 0 l - 1 f 1 ( X p [ i ] , U p [ i , · ] ) - x 0 ≤ 2 - p h l L F , R R . A24 Equation ( A23 ) also allows us to bound ‖ U p [ l , · ] - u 0 ‖ 1 . Since u 0 is piecewise smooth, ‖ U p [ l , · ] - u 0 ‖ 1 = 2 - p h Trap ( U p [ l , · ] - u 0 ) + O ( h ) . A25 Therefore, we need to bound the trapezoidal rule sum. We continue: 2 - p h Trap ( | U p [ l , · ] - u 0 | ) = 2 - p h Trap u 0 ( · ) + 2 - p h ∑ i = 1 l G ( · ) - u 0 ( · ) = 2 - p h Trap 2 - p h ∑ i = i 0 l - 1 G ( · ) = 2 - p h ∑ i = i 0 l - 1 2 - p h Trap | G ( · ) | , A26 where i 0 is 0 if the point index j ≥ l and i 0 = l - j if 0 ≤ j < l , and G ( · ) stands for f 2 ( · ) or 2 - p h Trap ( B ( · ) U p [ · ] ) – see Fig. 11 . Fig. 11. Open in a new tab An illustration to the proof of Lemma 5 , particularly to the calculation in ( A26 ). The large red squares mark the mesh with step h while the small black dots represent the mesh with step 2 - p h , where p = 3 . Here, l = 5 in ( A26 ). The large black dots mark the row with l = 5 . The summation in i of G ( · ) is done over the points of the same color arranged along the characteristics drawn with black lines. Al all colored dots lying on the y -axis, G ( · ) is 2 - p h Trap ( B ( · ) U p [ · ] ) , while at all other colored dots, G ( · ) is f 2 ( · ) Since ‖ f 2 ‖ 1 ≤ L F , R R , B ( · ) ≤ B max , and ‖ u 2 ‖ 1 ≤ R , we have ‖ U p [ l , · ] - u 0 ‖ 1 ≤ h 2 - p l max { B max , L F , R } R + O ( h ) ≤ A u h 2 - p l A27 for h ≤ 1 where A u : = max { B max , L F , R } R + 1 . Plugging ( A24 ) and ( A27 ) into ( A22 ) and computing the sum of the arithmetic series, we obtain the desired bound for the discrepancy in x between the coarse and fine mesh solutions over the first step in h : | X p [ 2 p ] - X 0 [ 1 ] | ≤ h 2 L F , R 2 p ( 2 p - 1 ) 2 - 2 p L F , R R + A u ≤ A 1 , x h 2 . A28 The discrepancy in u . Next, we need to show that U p [ 2 p , · ] - U 0 [ 1 , · ] 1 ≤ A 1 , u h 2 . A29 The values of U p and U 0 between the fine and coarse mesh points, respectively, are defined by linear interpolation. Therefore, the L 1 -norm is exactly equal to the quadrature by the trapezoidal rule on the fine mesh with step 2 - p h . We will the trapezoidal rule quadrature with step h and then estimate the correction due to U p . We start with the discrepancy at τ = 0 : U p [ 2 p , 0 ] - U 0 [ 1 , 0 ] = 2 - p h Trap B ( X p [ 2 p - 1 ] , · ) U p [ 2 p - 1 , · ] - h Trap B ( x 0 , · ) u 0 ( · ) ≤ 2 - p h Trap B ( X p [ 2 p - 1 ] , · ) U p [ 2 p - 1 , · ] - 2 - p h Trap B ( x 0 , · ) u 0 ( · ) + 2 - p h Trap B ( x 0 , · ) u 0 ( · ) - h Trap B ( x 0 , · ) u 0 ( · ) ≤ 2 - p h Trap B ( X p [ 2 p - 1 ] , · ) U p [ 2 p - 1 , · ] - 2 - p h Trap B ( x 0 , · ) U p [ 2 p - 1 , · ] + 2 - p h Trap B ( x 0 , · ) U p [ 2 p - 1 , · ] - 2 - p h Trap B ( x 0 , · ) u 0 ( · ) + 2 - p h Trap B ( x 0 , · ) u 0 ( · ) - h Trap B ( x 0 , · ) u 0 ( · ) To continue, we take into account the following facts: B is Lipschitz in x with constant L B and | X p [ 2 p - 1 ] - x 0 | ≤ h L F , R R according to ( A24 ); B ≤ B max and ‖ U p [ 2 p - 1 , · ] - u 0 ‖ 1 ≤ A u h according to ( A27 ); since B and u 0 are piecewise smooth, the difference between the trapezoidal rule quadrature with steps h and 2 - p h is O ( h ). Therefore, we conclude that if h is small enough, U p [ 2 p , 0 ] - U 0 [ 1 , 0 ] ≤ L B L F , R R h ‖ U p [ 2 p - 1 , · ] ‖ 1 + B max A u h + O ( h ) ≤ A 1 , 0 h . A30 Next, we calculate the discrepancies at a coarse mesh point with index [1, k ],  with k > 0 . This corresponds to the fine mesh index [ 2 p , 2 p k ] : U p [ 2 p , 2 p k ] - U 0 [ 1 , k ] = u 0 + 2 - p h ∑ i = 0 l - 1 f 2 ( · ) - u 0 - h f 2 ( x 0 , u 0 ( h ( k - 1 ) ) ) ≤ 2 - p h ∑ i = 0 2 p - 1 | f 2 [ p , i , 2 p ( k - 1 ) + i ] - f 2 ( x 0 , h ( k - 1 ) , u 0 ( h ( k - 1 ) ) | where we used a short-cut notation f 2 [ p , i , 2 p ( k - 1 ) + i ] : = f 2 ( X p [ i ] , 2 - p h ( 2 p ( k - 1 ) + i ) , U p [ i , 2 p ( k - 1 ) + i ] ) . A31 Therefore, h Trap ( | U p [ 2 p , · ] - U 0 [ 1 , · ] | ) ≤ h 2 A 1 , 0 h + h ∑ k ≥ 1 ′ U p [ 2 p , 2 p k ] - U 0 [ 1 , k ] ≤ h 2 2 A 1 , 0 + 2 - p h ∑ i = 0 2 p - 1 L F , R ‖ ( X p [ 2 p ] , U p [ 2 p , · ] ) - ( x 0 , u 0 ) ‖ X + O ( h ) ≤ A 1 ′ h 2 , where the sum with a prime above it indicates that the last term in it is divided by 2. This trapezoidal rule is exact for U 0 but is not exact for U p . However, on all intervals [ k h , ( k + 1 ) h ] with k ≥ 1 , U p [ 2 p , · ] can be extended to a piecewise smooth function thanks to such property of u 0 and f 2 . Therefore, the trapezoidal rule error on in most of these intervals will be O ( h 3 ) and in a few corresponding to jump discontinuities, O ( h 2 ) . The error on the first interval, [0, h ], will be O ( h 2 ) , because there U p can be extended continuously to [0, h ), but may have a discontinuity at h . Therefore, the overall error of the trapezoidal rule will be O ( h 2 ) . Therefore, we conclude that ( A29 ) holds. Moreover, U p [ 2 p , · ] can be extended into a piecewise smooth function. Step 2. The accumulation of the discrepancy between the solutions of coarse and fine meshes over a finite time interval. We will introduce the following notations e k , x : = X p [ 2 p k ] - X 0 [ k ] , e k , u [ · ] : = U p [ 2 p k , · ] - U 0 [ k , · ] . A32 We also will use a short-cut notation f 1 [ p , 2 p k ] : = f 1 ( X p [ 2 p k ] , U p [ 2 p k , · ] ) . A33 First, we compute e k + 1 , x given e k , x : e k + 1 , x : = X p [ 2 p ( k + 1 ) ] - X 0 [ k + 1 ] = X p [ 2 p k ] - X 0 [ k ] + 2 - p h ∑ l = 0 2 p - 1 f 1 [ p , 2 p k + l ] - f 1 [ 0 , k ] = e k , x + 2 - p h ∑ l = 0 2 p - 1 f 1 [ p , 2 p k + l ] - f 1 [ p , 2 p k ] + f 1 [ p , 2 p k ] - f 1 [ 0 , k ] . Using the Lipschitz continuity of f 1 and ( A21 ) we obtain | e k + 1 , x | ≤ | e k , x | + h L F , R | e k , x | + A 1 h 2 . A34 Second, we compute e k + 1 , u ( 0 ) given e k , u ( · ) : e k + 1 , u [ 0 ] : = 2 - p h Trap ( B [ p , 2 p k ] U p [ 2 p k , · ] ) - h Trap ( B [ 0 , 0 k ] U 0 [ k , · ] ) = 2 - p h Trap ( B [ p , 2 p k ] U p [ 2 p k , · ] ) - 2 - p h Trap ( B [ 0 , k ] U p [ 2 p k , · ] ) + 2 - p h Trap ( B [ 0 , k ] U p [ 2 p k , · ] ) - 2 - p h Trap ( B [ 0 , k ] U 0 [ k , · ] ) + 2 - p h Trap ( B [ 0 , k ] U 0 [ k , · ] ) - h Trap ( B [ 0 , k ] U 0 [ k , · ] ) , where B [ p , 2 p k ] : = B ( X p [ 2 p k ] , U [ 2 p k , · ] ) . Taking into account the facts listed in the bullet list above ( A30 ) we obtain | e k + 1 , u [ 0 ] | ≤ L B | e k , x | + ‖ e k , u ‖ 1 + B max ‖ e k , u ‖ 1 + O ( h ) . A35 Third, we compute e k + 1 , u ( j ) for j ≤ 1 referring to the coarse mesh in age. We will use a short-cut notation f 2 [ p , 2 p k , 2 p j ] : = f 2 ( X p [ 2 p k ] , 2 - p h 2 p k , U p [ 2 p k , 2 p j ] ) . A36 Thus, we calculate: e k + 1 , u [ j ] : = U p [ 2 p ( k + 1 ) , 2 p j ] - U 0 [ k + 1 , j ] = U p [ 2 p k , 2 p ( j - 1 ) ] + 2 - p h ∑ i = 0 2 p - 1 f 2 [ p , 2 p k + i , 2 p ( j - 1 ) + i ] - U 0 [ k , j - 1 ] - h f 2 [ 0 , k , j - 1 ] = e k , u [ j - 1 ] + 2 - p h ∑ i = 0 2 p - 1 f 2 [ p , 2 p k + i , 2 p ( j - 1 ) + i ] - f 2 [ p , 2 p k , 2 p ( j - 1 ) ] + 2 - p h ∑ i = 0 2 p - 1 f 2 [ p , 2 p k , 2 p ( j - 1 ) ] - f 2 [ 0 , k , j - 1 ] . Summing e k , u over the age range using the trapezoidal rule and taking into account that the numerical solutions can be extended to piecewise smooth functions, we get ‖ e k + 1 , u ‖ 1 ≤ h 2 L B ‖ e k ‖ X + h 2 B max + 1 ‖ e k , u ‖ 1 + h L F , R ‖ e k ‖ X + O ( h 2 ) , A37 where ‖ e k ‖ X : = | e k , x | + ‖ e k , u ‖ 1 . Adding ( A34 ) and ( A37 ) we obtain ‖ e k + 1 ‖ X ≤ ‖ e k ‖ X ( 1 + C 1 h ) + C 2 h 2 , A38 where C 1 and C 2 are appropriate positive constants. This implies that ‖ e k ‖ X ≤ ( 1 + C 1 h ) k ‖ e 0 ‖ X + ( 1 + C 1 h ) k - 1 1 + C 1 h - 1 C 2 h 2 = e C 1 h k ‖ e 0 ‖ X + e C 1 h k - 1 C 2 C 1 h . A39 Since h k ≤ h N 1 = T 1 and ‖ e 0 ‖ X = 0 , we obtain the desired bound: ‖ e k ‖ X ≤ e C 1 T 1 - 1 C 2 C 1 h . A40 □ Step 3. Taking the limit. The bound ( A40 ) shows that the sequence ( X p , U p ) of numerical solutions at any fixed time t < T 1 with steps 2 - p h is Cauchy in the Banach space X and hence converges to an element ( x , u ) in X . The x -component and u along the characteristics are continuously differentiable as immediately follows from time stepping. The function u ( t , τ ) may have discontinuities propagating characteristics emanating from the origin of the ( t , τ ) -plane and from the finite number of discontinuities of the initial age density u 0 – see the remark in Section 3.3 . Uniqueness Lemma 6 Under conditions of Theorem 1 , the solution to (7a)–(7h) is unique. Proof Let ( x , u ) and ( x ~ , u ~ ) be two solutions of (7a)–(7h) on the time interval [ 0 , T 1 ] defined in Section A.5 . We define a vector y ( t ) as follows: y ( t ) : = x ( t ) - x ~ ( t ) u ( t , 0 ) - u ~ ( t , 0 ) u ( t , t - t 0 ) - u ~ ( t , t - t 0 ) u ( t , t + τ 0 ) - u ~ ( t , t + τ 0 ) , where t 0 > 0 is the birth time, τ 0 > 0 is the initial age. We separate the age-zero term because it generates the initial conditions for all characteristics with t 0 ≥ 0 emanating from the t -axis. We introduce the norm ‖ y ( t ) ‖ Y = | x ( t ) - x ~ ( t ) | + | u ( t , 0 ) - u ~ ( t , 0 ) | + ∫ 0 ∞ | u ( t , τ ) - u ~ ( t , τ ) | d τ . A41 and a space Y consisting of all elements ( z , v 1 ( 0 ) , v 1 ( · ) , v 2 ( · ) ) with finite Y -norm: | z | + | v 1 ( 0 ) | + ‖ v 1 ‖ 1 + ‖ v 2 ‖ 1 < ∞ . The Y -space is Banach as it is a direct product of a Banach space and an intersection of Banach spaces. As in the X -space, we consider a ball of a large radius R in the Y -space and assume that the solutions ( x , u ) and ( x ~ , u ~ ) lie in this ball. Next, we consider a Fréchet derivative of y ( t ): y ′ ( t ) = f 1 ( x ( t ) , u ( t , · ) ) - f 1 ( x ~ ( t ) , u ~ ( t , · ) ) ∫ 0 ∞ B ( x ( t ) , τ ) u t ( t , τ ) d τ - ∫ 0 ∞ B ( x ~ ( t ) , τ ) u ~ t ( t , τ ) d τ f 2 ( x ( t ) , t - t 0 , u ( t , t - t 0 ) ) - f 2 ( x ~ , t - t 0 , ( t ) , u ~ ( t , t - t 0 ) ) f 2 ( x ( t ) , t + τ 0 , u ( t , t - t 0 ) ) - f 2 ( x ~ , t + τ 0 , ( t ) , u ~ ( t , t + τ 0 ) ) . A42 Using the PDE for the predator age density (5b), we rewrite integrals in the second component of y ′ ( t ) as follows: ∫ 0 ∞ B ( x ( t ) , τ ) u t ( t , τ ) d τ = - ∫ 0 ∞ B ( x ( t ) , τ ) u τ ( t , τ ) d τ - ∫ 0 ∞ B ( x ( t ) , τ ) μ ( x ( t ) , τ ) u ( t , τ ) d τ = - B ( x ( t ) , 0 ) u ( t , 0 ) + ∫ 0 ∞ B τ ( x ( t ) , τ ) u ( t , τ ) d τ - ∫ 0 ∞ B ( x ( t ) , τ ) μ ( x ( t ) , τ ) u ( t , τ ) d τ . Taking into account that u ( t , τ ) has a compact support in τ at any t < ∞ as u 0 ( τ ) is compactly supported, and that ( x , u ) and ( x ~ , u ~ ) lie in the ball of radius R in the Y -space, we obtain the following bound for the second component of y ′ ( t ) ∫ 0 ∞ B ( x , τ ) u t ( t , τ ) d τ - ∫ 0 ∞ B ( x ~ ( t ) , τ ) u ~ t ( t , τ ) d τ ≤ | B ( x , 0 ) | | u ( t , 0 ) - u ~ ( t , 0 ) | + | B ( x , 0 ) - B ( x ~ , 0 ) | | u ~ ( t , 0 ) | + ∫ 0 ∞ | B τ ( x , τ ) | | u ( t , τ ) - u ~ ( t , τ ) | d τ + ∫ 0 ∞ | B τ ( x , τ ) - B τ ( x ~ , τ ) | | u ~ ( t , τ ) | d τ + ∫ 0 ∞ | B ( x , τ ) μ ( x , τ ) | | u ( t , τ ) - u ~ ( t , τ ) | d τ + ∫ 0 ∞ | B ( x , τ ) μ ( x , τ ) - B ( x ~ , τ ) μ ( x ~ , τ ) | | u ~ ( t , τ ) | d τ ≤ | B ( x , 0 ) | | u ( t , 0 ) - u ~ ( t , 0 ) | + L B τ | x - x ~ | R + B max μ max ‖ u ( t , · ) - u ~ ( t , · ) ‖ 1 + L B μ | x - x ~ | R , where L B τ and L B μ are the Lipschitz constants for B τ and B μ respectively. Note that B τ ( x , τ ) is Lipschitz in x because the smoothed indicator function has finite derivative at all τ . Taking Lipschitz continuity of f 1 and f 2 into account and using the bound for the second component of y ′ ( t ) just derived, we claim that ‖ y ′ ( t ) ‖ Y ≤ L ‖ y ( t ) ‖ Y A43 for an appropriate Lipschitz constant L . Next, we fix a small h > 0 . By the triangle inequality, ‖ y ( t + h ) ‖ Y ≤ ‖ y ( t + h ) - y ( t ) ‖ Y + ‖ y ( t ) ‖ Y . Hence, ‖ y ( t + h ) ‖ Y - ‖ y ( t ) ‖ Y h ≤ ‖ y ( t + h ) - y ( t ) ‖ Y h . The same is true if h < 0 and t > - h . Together with ( A43 ) this implies that d dt ‖ y ‖ Y ≤ ‖ y ′ ( t ) ‖ ≤ L ‖ y ‖ Y . A44 Therefore, by Grönwall’s Inequality, 0 ≤ ‖ y ( t ) ‖ Y ≤ ‖ y ( 0 ) ‖ e Lt . A45 Since y ( 0 ) = 0 , it remains zero for all times. Hence the solution to (7a)–(7h) is unique. □ Positivity Positivity of x ( t ). Lemma 7 Under the conditions of Theorem 1 , the x -component of the solution is positive provided that x (0) is positive. Proof Equation (7a) for x can be rewritten as follows: x ′ ( t ) = x ( t ) f ( x ( t ) , u ( t , · ) ) , A46 where f ( t ) = r - a x ( t ) + s ∫ 0 ∗ u ( t , τ ) d τ - b ∫ τ ∗ τ max ( T ) u ( t , τ ) d τ . A47 Therefore, we can introduce z ( t ) : = log x ( t ) and write the following ODE for z : z ′ = f ( e z , u ( t , · ) ) . A48 If we replace x with z in (7a)–(7h), we can discretize it and conduct the existence and uniqueness proofs essentially repeating the arguments in the previous sections of Appendix A. Importantly, the initial condition for z , z ( 0 ) = log x ( 0 ) is defined. Hence, a solution ( z , u ) exists. Then x ( t ) = e z ( t ) also exists and is positive. □ Positivity of u ( t , τ ) . Lemma 8 Let u 0 ( τ ) be positive on [ 0 , τ 0 , max ) . Under the conditions of Theorem 1 , the u -component of the solution is positive in the trapezoidal ( t , τ ) -region T : = { ( t , τ ) | 0 ≤ t ≤ T , 0 ≤ τ < τ 0 , max + t } . A49 Proof One can argue that u ( t , t + τ 0 ) > 0 and u ( t , t - t 0 ) > 0 as soon as u 0 ( τ 0 ) > 0 and u ( t 0 , 0 ) > 0 in the same manner as in the proof of Lemma 7 . It is given that u 0 ( τ ) > 0 on [ 0 , τ 0 , max ) . Therefore, the only opportunity for u to become nonpositive in T is to acquire a nonpositive initial condition at ( t 0 , 0 ) for some t 0 ≥ 0 . This means that the birth integral must become nonpositive. Let t 0 ≥ 0 be the smallest time at which u ( t 0 , 0 ) = ∫ 0 ∞ B ( x ( t 0 ) , τ ) u ( t 0 , τ ) d τ ≤ 0 . This mean that u ( t 0 , τ ) ≤ 0 at some set of τ of positive measure inside [ 0 , τ 0 , max + t 0 ) . Let τ 1 ∈ [ 0 , τ 0 , max + t 0 ) be such that u ( t 0 , τ 1 ) ≤ 0 . Then, since u preserves its sign along the characteristics, u must be negative along the whole characteristic passing through ( t 0 , τ 1 ) . If this characteristic hits the τ -axis first, we come to a contradiction with the assumption that u 0 > 0 on [ 0 , τ 0 , max ) . If this characteristic hits the t -axis first, we arrive at a contradiction with the assumption that t 0 is the earliest moment of time when u ( t , τ ) became nonpositive on τ ∈ [ 0 , τ 0 , max + t ) . Hence, u ( t , τ ) is positive in T . □ Linear Discriminant Analysis Let X be a n × d matrix whose rows are d -dimensional data points. The dataset X consists of c > 1 categories. Let I i denote the set of indices of data from category i , i = 1 , … , c , and | I i | = n i , I i ∩ I j = ∅ , I 1 ∪ … ∪ I c = { 1 , … , n } . We seek a d × d 1 matrix W , d 1 ≤ c - 1 , that maps the data onto a d 1 -dimensional space: Y = X W , or , y k = W ⊤ x k , k = 1 , … n , B50 such that images of the data from different categories under this mapping are separated as much as possible while data from the same categories are clustered as much as possible. In (B50), x i ∈ R d is the k th row of X written as a column vector. Likewise is y k . The means for each category in the spaces R d and R d 1 are, respectively, m i = 1 n i ∑ k ∈ I i x k , and m ~ i = 1 n i ∑ k ∈ I i y k = W ⊤ m i , i = 1 , … , c . B51 To describe the data variation within and between the categories, the within-class and between-class scatter matrices S w and S b are introduced. These matrices are of size d × d . The within-class scatter matrix is defined as S w : = ∑ i = 1 c S i , B52 where S i is the scatter matrix of category i defined as S i : = ∑ k ∈ I i ( x k - m i ) ( x k - m i ) ⊤ ≡ X I i , : - 1 n i × 1 m i ⊤ ⊤ X I i , : - 1 n i × 1 m i ⊤ . B53 The between-class scatter matrix is defined as S b : = ∑ i = 1 c n i ( m i - m ) ( m i - m ) ⊤ , where m : = 1 n ∑ i = 1 c n i m i ≡ 1 n ∑ k = 1 n x k B54 is the overall mean. The rank of S b is at most c - 1 . A similar notation with tilde on top is used for the mapped data. The within- and between-class scatter matrices for the mapped data are: S ~ w = W ⊤ S w W , S ~ b = W ⊤ S b W . B55 The task of LDA is to find a d × d 1 matrix W , d 1 ≤ c - 1 , mapping the data within each category into clusters and mapping the clusters corresponding to different categories as far as possible from each other. Respectively, the LDA objective function is defined as J ( w ) = w ⊤ S b w w ⊤ S w w . B56 The function J ( w ) is maximized by solving the generalized eigenvalue problem. Indeed, the gradient of J ( w ) is zero if and only if S b w = J ( w ) S w w . B57 Therefore, to find the desired projection one needs to solve (B57) and compose the projection matrix W out of d 1 eigenvectors corresponding to d 1 largest eigenvalues. Author Contributions Luis Suarez: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing – original draft, Visualization. Maria Cameron: Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing – review&editing, Visualization, Supervision, Funding Acquisition. William Fagan: Conceptualization, Writing – review&editing, Supervision. Doron Levy: Conceptualization, Writing – review& editing, Supervision, Project administration, Funding Acquisition. Funding AFOSR MURI grant FA9550-20-1-0397 (MC); Simons Foundation Award 848629 (DL); U.S. National Science Foundation Award DMS2451241 (WFF). Data Availability Due to the volume of the simulation data, it will be made available upon request. Declarations Conflict of interest/Competing interests The authors have no conflict of interest. Ethics approval and consent to participate The authors comply with ethics rules. Consent for publication All authors consent for publication. Materials availability Not applicable. Code availability Codes are available on GitHub (Suarez and Cameron 2025 ). Footnotes Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. L.C. Suarez, M.K. Cameron:These authors contributed equally to this work. References Audze P, Eglãjs V (1977) New approach to the design of multifactor experiments. Problems of Dynamics and Strengths 35:104–107 (in Russian) Adams RA, Fournier JJF (2003) Sobolev Spaces. Elsevier, Amsterdam, London Allee WC, Park O, Emerson AE, Schmidt TPKP (1949) Principles of Animal Ecology. Saunders Company, Philadelphia 10.5962/bhl.title.7325 Agutter PS, Wheatley DN (2008) Thinking About Life: The History and Philosophy of Biology and Other Sciences. Springer, Dordrecht, London. https://archive.org/details/thinkingaboutlif0000agut Cantrell RS, Cosner C, Fagan WF (2012) The implications of model formulation when transitioning from spatial to landscape ecology. Math Biosci Eng 9(1):27–60. 10.3934/mbe.2012.9.27 [ DOI ] [ PubMed ] [ Google Scholar ] Cushing JM, Saleem M (1982) A predator prey model with age structure. J Math Biol 14:231–250. 10.1007/BF01832847 [ DOI ] [ PubMed ] [ Google Scholar ] Cushing JM (1984) Existence and stability of equilibria in age-structured population dynamics. J Math Biol 20:259–276. 10.1007/BF00275988 [ Google Scholar ] Cushing JM (1986) Periodic McKendrick equations for age-structured population growth. In: WITTEN, M. (ed.) Hyperbolic Partial Differential Equations, pp. 513–526. Pergamon Press, Oxford, UK. 10.1016/B978-0-08-034313-6.50013-4 . https://www.sciencedirect.com/science/article/pii/B9780080343136500134 Cushing JM (1992) A size-structured model for cannibalism. Theor Popul Biol 42(3):347–361. 10.1016/0040-5809(92)90020-T [ Google Scholar ] Duda RO, Hart PE, Stork DG (2001) Pattern Classification, 2nd edn. John Wiley and Sons Inc, New York [ Google Scholar ] Elton CS (1927) Animal Ecology. Text-books of animal biology, ed, by Julian S. Huxley. Sidgwick & Jackson, London, UK. http://books.google.com/books?id=14BOAQAAIAAJ Fouilloux CA, Fromhage L, Valkonen JK, Rojas B (2022) Size-dependent aggression towards kin in a cannibalistic species. Behav Ecol 33(3):582–591. 10.1093/beheco/arac020 [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] Fagan WF, Odell GM (1996) Size-dependent cannibalism in praying mantids: Using biomass flux to model size-structured populations. Am Nat 147(2):230–268. 10.1086/285848 [ Google Scholar ] Iannelli M (1995) Mathematical Theory of Age-structured Population Dynamics. Giardini editori e stampatori, Pisa https://api.semanticscholar.org/CorpusID:117215818 Iman RL, Helton JC, Campbell JE (1981) An Approach to Sensitivity Analysis of Computer Models: Part I—Introduction, Input Variable Selection and Preliminary Variable Assessment. J Qual Technol 13(3):174–183. 10.1080/00224065.1981.11978748 [ Google Scholar ] Köster FW, Möllmann C (2000) Trophodynamic control by clupeid predators on recruitment success in baltic cod? ICES J Mar Sci 57:310–323. 10.1006/jmsc.1999.0528 Lehtinen SO (2021) Ecological and evolutionary consequences of predator-prey role reversal: Allee effect and catastrophic predator extinction. Journal of Theoretical Biology 510, 110542 10.1016/j.jtbi.2020.110542 Li J, Liu X, Wei C (2022) The impact of role reversal on the dynamics of predator-prey model with stage structure. Appl Math Model 104:339–357. 10.1016/j.apm.2021.11.029 [ Google Scholar ] Lotka AJ (1925) Elements of Physical Biology. Nature 116:461. 10.1038/116461b0 [ Google Scholar ] McKay MD, Beckman RJ, Conover WJ (1979) A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21(2):239–245. 10.2307/1268522 [ Google Scholar ] Mohr M, Barbarossa MV, Kuttler C (2014) Predator-prey interactions, age structures and delay equations. Math Model Nat Phenom 9(1):92–107. 10.1051/mmnp/20149107 [ Google Scholar ] McKendrick AG, Pai MK (1912) XLV.—The Rate of Multiplication of Micro-organisms: A Mathematical Study. https://api.semanticscholar.org/CorpusID:86967998 Mishra P, Ponosov A, Wyller J (2024) On the dynamics of predator-prey models with role reversal. Physica D 461:134100. 10.1016/j.physd.2024.134100 [ Google Scholar ] Nakazawa T (2015) Ontogenetic niche shifts matter in community ecology: a review and future perspectives. Popul Ecol 57(2):347–354. 10.1007/s10144-014-0448-z [ Google Scholar ] Neuenfeldt S, Köster FW (2000) Trophodynamic control on recruitment success in Baltic cod: the influence of cannibalism. J Mater Sci 57:300–309. 10.1006/jmsc.1999.0528 [ Google Scholar ] Nocedal J, Wright SJ (2006) Numerical Optimization, 2nd edn. Springer, New York, NY, USA [ Google Scholar ] Perthame B (2006) Transport Equations in Biology. Frontiers in Mathematics. Birkhäuser Basel. https://books.google.com/books?id=1vZlBfhCEgMC Pavel NH, Motreanu D (1999) Tangency, Flow Invariance for Differential Equations, and Optimization Problems. CRC Press, New York, NY, USA [ Google Scholar ] Polis GA (1988) Exploitation competition and the evolution of interference, cannibalism, and intraguild predation in age/size-structured populations. In: Size Structured Populations: Ecology and Evolution. https://api.semanticscholar.org/CorpusID:81254128 Polis GA (1991) Complex trophic interactions in deserts: An empirical critique of food-web theory. Am Nat 138(1):123–155. 10.1086/285208 [ Google Scholar ] Saitā, Y (1986) Prey kills predator: Counter-attack success of a spider mite against its specific phytoseiid predator. Experimental & Applied Acarology 2, 47–62 10.1007/BF01193354 Suarez L, Cameron M (2025) GitHub: Predator-Prey. A predator-prey model with an age-structured role reversal. https://github.com/mar1akc/Predator-Prey/tree/main Volterra L (1926) Variazioni e Fluttuazioni del Numero d’Individui in Specie Animali Conviventi. Città di Castello: Società anonima tipografica “Leonardo da Vinci ” https://liberliber.it/autori/autori-v/vito-volterra/variazioni-e-fluttuazioni-del-numero-dindividui-in-specie-animali-conviventi/ Werner EE, Gilliam JF (1984) The ontogenetic niche and species interactions in size-structured populations. Annu Rev Ecol Evol Syst 15:393–425. 10.1146/annurev.es.15.110184.002141 [ Google Scholar ] Li Q, Lin B, Ren W (2019) Computing committor functions for the study of rare events using deep learning. J Chem Phys 151:054112. 10.1063/1.5110439 [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Data Availability Statement Due to the volume of the simulation data, it will be made available upon request. Articles from Journal of Mathematical Biology are provided here courtesy of Springer ACTIONS View on publisher site PDF (2.6 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top