Research article

Global dynamics of a discrete two-population model for flour beetle growth


  • Received: 29 April 2025 Revised: 27 May 2025 Accepted: 09 June 2025 Published: 25 June 2025
  • The cannibalistic behavior of Tribolium has been extensively researched, revealing instances of chaotic dynamics in laboratory environments for Tribolium castaneum. The well-established Larvae-Pupae-Adult (LPA) model has been instrumental in understanding the conditions that lead to chaos in flour beetles (genus: Tribolium). In response to new experimental observations showing a decline in the pupae population in Tribolium confusum, we proposed and analyzed a simplified two-stage Larvae-Adult (LA) model. This model integrated the pupae population within the larval group, similar to that of the original LPA model, with development transitions governed by internal rates. By applying the model to time-series data, we demonstrated its effectiveness in capturing short-term population fluctuations in T. confusum. We established the model's positivity and boundedness, perform stability analyses of both trivial and positive steady states, and explored bifurcations and steady-state behavior through numerical simulations. We proved global stability for the extinction and positive steady states and observed additional restrictions required for stability compared to the LPA model. Our results indicated that while chaos was a possible outcome, it was infrequent within the practical parameter ranges observed, with environmental changes related to media and nutrient alterations being more likely triggers.

    Citation: Samantha J. Brozak, Kamrun N. Keya, Denise Dengi, Sophia Peralta, John D. Nagy, Yang Kuang. Global dynamics of a discrete two-population model for flour beetle growth[J]. Mathematical Biosciences and Engineering, 2025, 22(8): 1980-1998. doi: 10.3934/mbe.2025072

    Related Papers:

    [1] Horst R. Thieme . Discrete-time population dynamics on the state space of measures. Mathematical Biosciences and Engineering, 2020, 17(2): 1168-1217. doi: 10.3934/mbe.2020061
    [2] Francesca Verrilli, Hamed Kebriaei, Luigi Glielmo, Martin Corless, Carmen Del Vecchio . Effects of selection and mutation on epidemiology of X-linked genetic diseases. Mathematical Biosciences and Engineering, 2017, 14(3): 755-775. doi: 10.3934/mbe.2017042
    [3] Yicang Zhou, Zhien Ma . Global stability of a class of discrete age-structured SIS models with immigration. Mathematical Biosciences and Engineering, 2009, 6(2): 409-425. doi: 10.3934/mbe.2009.6.409
    [4] Qianhong Zhang, Fubiao Lin, Xiaoying Zhong . On discrete time Beverton-Holt population model with fuzzy environment. Mathematical Biosciences and Engineering, 2019, 16(3): 1471-1488. doi: 10.3934/mbe.2019071
    [5] Yanjie Hong, Wanbiao Ma . Sufficient and necessary conditions for global attractivity and stability of a class of discrete Hopfield-type neural networks with time delays. Mathematical Biosciences and Engineering, 2019, 16(5): 4936-4946. doi: 10.3934/mbe.2019249
    [6] Yufei Wang, Huidong Cheng, Qingjian Li . Dynamic analysis of wild and sterile mosquito release model with Poincaré map. Mathematical Biosciences and Engineering, 2019, 16(6): 7688-7706. doi: 10.3934/mbe.2019385
    [7] Chang-Yuan Cheng, Shyan-Shiou Chen, Xingfu Zou . On an age structured population model with density-dependent dispersals between two patches. Mathematical Biosciences and Engineering, 2019, 16(5): 4976-4998. doi: 10.3934/mbe.2019251
    [8] Shuixian Yan, Sanling Yuan . Critical value in a SIR network model with heterogeneous infectiousness and susceptibility. Mathematical Biosciences and Engineering, 2020, 17(5): 5802-5811. doi: 10.3934/mbe.2020310
    [9] Gheorghe Craciun, Stefan Muller, Casian Pantea, Polly Y. Yu . A generalization of Birchs theorem and vertex-balanced steady states for generalized mass-action systems. Mathematical Biosciences and Engineering, 2019, 16(6): 8243-8267. doi: 10.3934/mbe.2019417
    [10] Ning Yu, Xue Zhang . Discrete stage-structured tick population dynamical system with diapause and control. Mathematical Biosciences and Engineering, 2022, 19(12): 12981-13006. doi: 10.3934/mbe.2022606
  • The cannibalistic behavior of Tribolium has been extensively researched, revealing instances of chaotic dynamics in laboratory environments for Tribolium castaneum. The well-established Larvae-Pupae-Adult (LPA) model has been instrumental in understanding the conditions that lead to chaos in flour beetles (genus: Tribolium). In response to new experimental observations showing a decline in the pupae population in Tribolium confusum, we proposed and analyzed a simplified two-stage Larvae-Adult (LA) model. This model integrated the pupae population within the larval group, similar to that of the original LPA model, with development transitions governed by internal rates. By applying the model to time-series data, we demonstrated its effectiveness in capturing short-term population fluctuations in T. confusum. We established the model's positivity and boundedness, perform stability analyses of both trivial and positive steady states, and explored bifurcations and steady-state behavior through numerical simulations. We proved global stability for the extinction and positive steady states and observed additional restrictions required for stability compared to the LPA model. Our results indicated that while chaos was a possible outcome, it was infrequent within the practical parameter ranges observed, with environmental changes related to media and nutrient alterations being more likely triggers.



    Tribolium have infested stored grains for thousands of years, making them an ideal model organism for study due to their proximity to humans and adaptation to the experimental environment [1]. Also known as flour beetles, they have contributed extensively to scientific advancements in evolution, reproduction, and, perhaps surprising to some, mathematics [1,2,3]. An initial claim to fame was led by Thomas Park and colleagues, whose two-species competition experiments with seemingly random outcomes spurred mathematical interest [4,5,6,7]. A second wave of fame, this time led by the Beetle Team, used flour beetles in the search for chaos in a laboratory population [2,3,8,9]. The team successfully induced chaos in Tribolium castanaeum in their landmark studies in the late 1990s and inspired recent work by [10].

    Aside from their adaptation to the laboratory, flour beetle biology contains many nonlinearities that contribute to their scientific interest. Tribolium are holometabolous and thus go through four major life stages: Egg, larva, pupa, and adult [1]. Egg to adult development takes approximately 30 days in optimal conditions (around 30 C) [1]. Eggs develop over the course of approximately three days, with twenty days spent as larvae and four days spent as pupae [1]. After pupation is the callow stage, in which the new adult beetle is still reproductively immature and does not have a hardened exoskeleton for three to four days. Sclerotized adults, or adults with a hardened exoskeleton, are dark brown in color and 4 mm in length [1]. Adults have the ability to fly but do not tend to do so in laboratory settings [1]. The relatively quick reproduction time and ease of handling makes Tribolium an excellent model organism for scientific study.

    Flour beetles regulate their population via cannibalism, with larvae mostly consuming eggs and adults consuming eggs and pupae [1,11,12]. Cannibalism may provide an additional way of obtaining nutrients in a nutrient-poor environment [13], and thus facilitate colonization of a new environment [14]. It is thought that cannibalism drives the oscillations observed in flour beetles [1]. These oscillations captured researchers' interest when Robert May spurred a hunt for chaos in nature [8] and motivated researchers' choices to use flour beetles [2,3].

    The LPA model is described by the following discrete-time system [15]:

    L(t+1)=bA(t)ec1L(t)c2A(t),P(t+1)=(1μl)L(t),A(t+1)=P(t)ec3A(t)+(1μa)A(t). (1.1)

    The number of eggs that hatch into larvae in the absence of cannibalism is given by b>0. Parameters μl,μa(0,1) denote losses due natural mortality (excluding cannibalism) for larvae and adults, respectively. Cannibalism is modeled using a binomial distribution and approximated using the exponential where ci, i=1,2,3 represent the intensity of cannibalism. The LPA model (1.1) has been well-characterized analytically, from local and global stability to cycles and cycle chains [16,17,18,19,20,21,22]. Furthermore, the model has been validated rigorously using a variety of statistical methodologies [9,23,24].

    New experimental data has inspired more work with Tribolium. Brozak et al. [10] used new T. confusum data to develop a model that stratified adults by reproductive capacity, called the LPAA model. The authors found that the LPA model (1.1) could not adequately describe their data, likely due to differences in experimental and censusing protocols, including species used. The results suggested that chaos was not an inherent property of these beetle populations and must be induced [10].

    Using new, longer-term experimental data for T. confusum, we develop and parameterize a two-stage discrete-time model, which we refer to as the LA model. We derive thresholds for local and global stability of the extinction and positive steady states, and find that the conditions for global stability are stricter than previous models.

    We propose the following model, which we refer to as the LA model:

    L(t+1)=bA(t)ec1A(t)+(1μl)L(t), (2.1)
    A(t+1)=μlL(t)ec2A(t)+(1μa)A(t), (2.2)

    where L(t) describes the total population of larvae and pupae at time t, and A(t) is the number of adults at time t. Similar to previous work, one timestep corresponds to two weeks [10,15]. Parameter b>0 denotes larval recruitment; c1,c2(0,1) are the cannibalism intensity from adults consuming eggs and adults consuming pupae, respectively. Like the LPA model, we assume that the contacts between adults and an egg or pupa occur at random. The probability that the egg or pupa survives is given by ec1A(t) and ec2A(t) [23]. We note that, while one could include the exponential cannibalism term in the losses for Eq (2.1), these can be indirectly accounted for in the recruitment term, especially so since the juvenile populations are low in our data. Here, μl(0,1) is the proportion of L(t) that transition to adulthood, while 1μl remain in L(t) to continue pupation. Parameter μa(0,1) is the mortality probability for adults.

    In the LA model, larvae and pupae are combined into a single population category due to the improved fit against the observed decline in pupae numbers in our experimental data (discussed in Section 3). This grouping of the juvenile stages is similar to that of the original LPA work in which large (non-feeding) larvae were grouped with pupae, although T. confusum larvae continue to feed in later instars [25]. The life cycle as included in the model is shown in Figure 1.

    Figure 1.  Life cycle of T. confusum in the absence of cannibalism. Arrows indicate the transitions between life stages. The LA model incorporates larvae and pupae as one stage (indicated by the box) as well as adults.

    We compare the LA model (2.1) and (2.2) to new experimental data looking at the effects of nutrient ratios on populations of T. confusum. The experimental protocol is similar to that of Brozak et al. [10] where insects were counted and reported every two weeks, except with higher ratios of nitrogen (N) and phosphorus (P). Additional details are available in Section ??? of this manuscript. The decline of all populations in the study is observed, and there are several hypotheses for why this may occur. Higher concentrations of nitrogen and phosphorus in the environment may have reached toxic levels, but the decline is also observed in the control group. A large larval die off was observed at week 16 potentially due to loss of power in the environmental chambers, but large adult die offs were observed and recorded prior to this event. Specifically, the decreasing number of pupae prompted the introduction of the LA model (2.1) and (2.2), instead of the LPA [26] and LPAA [10] models.

    The fitting process is carried out for all experimental groups, as shown in Table 1. The sum of larval and pupal data is used for parameterizing the L(t) while the total adult population (given by unsclerotized adults and sclerotized adults) is used for fitting A(t). Parameter set 1 is obtained from the aggregated data (without the control groups), where the data is combined using a pointwise mean. Parameter sets 2 and 3 are fit to the data from the cups with 5% and 10% of added nitrogen (N) by mass, respectively. Parameter set 4 is calculated from the data from cups with 5% of added nitrogen (N) and 5% of added phosphorus (P), and parameter set 5 is based on the carbon (C) control data (no added N or P). The naming convention for the groups will be the element followed by its ratio (for example, N5P5 indicates 5% nitrogen and 5% phosphorus).

    Table 1.  Best fit parameters of the model.
    Data type b μa c1 c2 μl L Error A Error
    1 Aggregated 52.7 0.1992 0.0540 3×104 0.3203 27% 10%
    2 N5 30.30 0.2574 0.0251 1.34×105 0.6044 18% 20%
    3 N10 27.7 0.3241 0.0377 2.18×104 0.5480 21% 19%
    4 N5P5 46 0.3749 0.0309 1.03×105 0.4566 20% 18%
    5 C 39.4 0.2297 0.0316 4.40×104 0.4347 18% 10%

     | Show Table
    DownLoad: CSV

    MATLAB's "fmincon" and "MultiStart" functions are used to find the best approximated parameters and fittings. Due to the low cannibalism rate, the the bounds for c1 and c2 are set at [0,0.2]. The following normalized relative error formula is used for optimization:

    t|L(t)Ldata(t)|t|L(t)+Ldata(t)|+t|A(t)Adata(t)|t|A(t)+Adata(t)|, (3.1)

    where, for X={L,A}, X(t) corresponds to the model output at time t, and Xdata(t) is the experimental data at time t. This is a different form than the typical relative error formula; this formulation ensures that lower magnitude data is not overweighted.

    Comparisons between the LPA, LPAA, and the LA model are shown in Figure 2, with the LA model performing the best on the time series. The LPA and LPAA models failes to capture the sustained oscillations in the larvae. Some best fit results are presented in Figures 3 and 4. For Figure 3 (left), the best approximations of parameters are shown in parameter set 1 in Table 1, and the best approximated parameters for Figure 3 (right) are presented in parameter set 5 of Table 1. In Figure 4, fitting results for data groups with added nutrients are presented, with all the respective parameter sets shown in Table 1. Each figure exhibits two distinct phases. The first is a transient phase occurring in the first 8–10 weeks, during which there was a notable peak in the larvae population followed by a sudden decline. The remaining data shows a similar structure to the transient phase, but displays a recurring pattern for the rest of the time period referred to as the intermediate phase.

    Figure 2.  Performance of the (a) LPA, (b) LPAA, and (c) LA model on aggregated data. The LA model demonstrates the closest fit to the real data, capturing the observed population dynamics.
    Figure 3.  Model fitting for aggregated data (left) and control group data (right).
    Figure 4.  Model fittings for data group N5 (left) and data group N5P5 (right) are shown. Best approximated parameters for each data set are provided in Table 1.

    In this section, we show that solution of the LA model does not extend beyond the region Ω. The following lemma defines Ω and provides a straightforward proof of the boundedness of the LA model. To demonstrate this boundedness, we utilize basic iteration methods.

    Lemma 4.1. The solution of system (2.1) and (2.2) is eventually uniformly bounded, and the following region is positively invariant with respect to system (2.1) and (2.2):

    Ω={(L,A)R2+:Lbec1μl,Abec1μa}.

    Proof. We first show that solutions are eventually uniformly bounded. This is equivalent to showing that the lim sup of L and A are bounded by constants independent of initial values.

    Assume that the initial population size of larvae and adult are L(0)=L00 and A(0)=A0>0. From system (2.1) and (2.2), we have,

    L(t)0andA(t)>0for allt0.

    Let,

    lt=max{L(i);i=0,1,...,t},at=max{A(i);i=0,1,...,t}.

    Here, lt and at are increasing functions of t. We also define m=max{bxec1x;x>0}=bec1. From (2.1), we see

    L(t+1)=bA(t)ec1A(t)+(1μl)L(t)m+(1μl)ltm+(1μl)lt+1.

    Therefore,

    L(i)m+(1μl)lim+(1μl)li+1,for all i=0,1,,t.

    Hence,

    li+1m+(1μl)li+1li+1mμl,

    which implies limtsupL(t+1)mμl=bec1μl.

    Since L(t) is bounded for all t0, using a similar argument in (2.2), we have

    A(t+1)=μlL(t)ec2A(t)+(1μa)A(t)mec2A(t)+(1μa)A(t)m+(1μa)A(t)m+(1μa)at.

    Henceforth, we have

    at+1m+(1μa)at+1at+1mμa.

    Therefore, limtsupA(t+1)mμa=bec1μa.

    As for the invariance of Ω, it is easy to see that it follows from a straightforward mathematical induction on li and ai.

    The LA model (2.1) and (2.2) admits two steady states. Suppose (ˉL,ˉA) is an arbitrary steady state of the LA model. Then L(t+1)=L(t)=ˉL and A(t+1)=A(t)=ˉA, which yields

    ˉL=bˉAec1ˉA+(1μl)ˉL,ˉA=μlˉLec2ˉA+(1μa)ˉA.

    Solving this system gives an extinction state, E0=(0,0), and a positive steady state E1=(L,A), where L and A satisfy

    L=bμlAec1A, (4.1)
    A=1c1+c2ln(bμa). (4.2)

    The steady state A exists and is unique when ln(bμa)>0, or bμa>1. Therefore, we can define R0=bμa>1. This quantity represents the average number of offspring from a single adult that reach adulthood.

    We analyze the local stability of the extinction and positive steady states of system (2.1) and (2.2) using the classic method of analyzing the Jacobian. From [27] we restate the following lemma:

    Lemma 4.2 (Edelstein-Keshet [27]). Suppose x and y are two species population variables. Consider the system of equations

    xn+1=f(xn,yn), (4.3)
    yn+1=g(xn,yn), (4.4)

    where f and g are nonlinear functions, and (ˉx,ˉy) is a steady state of (4.3) and (4.4). The Jacobian matrix of this system is given by

    J=[a11a12a21a22]=[fx|(ˉx,ˉy)fy|(ˉx,ˉy)gx|(ˉx,ˉy)gy|(ˉx,ˉy)].

    The characteristic polynomial is then λ2βλ+γ=0, where β=a11+a22 and γ=a11a22a12a21. If |β|<1+γ<2, then |λi|<1, i={1,2}, and so (ˉx,ˉy) is locally asymptotically stable.

    In the following theorem, we derive thresholds for local stability of the extinction and positive steady states.

    Theorem 4.3. The following properties hold for the LA model (2.1) and (2.2):

    1) Extinction steady state E0=(0,0) is locally asymptotically stable when R0<1 and unstable when R0>1.

    2) Positive steady state E1=(L,A) defined in (4.1) and (4.2) exists when R0>1.

    (a) If μl>c2c1+c2, then E is locally asymptotically stable for

    1<R0<exp(μl+μaμa(μlc2c1+c2)).

    (b) If μl<c2c1+c2, then E is locally asymptotically stable for R0>1.

    Proof. The Jacobian matrix of system (2.1) and (2.2) around the general steady state is

    JLA=[1μl(1c1ˉA)bec1ˉAμlec2ˉA1μac2μlˉLec2ˉA]. (4.5)

    For the extinction steady state E0=(0,0), we have

    JLA(0,0)=[1μlbμl1μa].

    For local stability, we need to satisfy |β|<1+γ<2 provided in Lemma 4.2. We have

    β=2μlμa0,γ=(1μl)(1μa)μlb.

    We have

    1+γ=β+μl(μab)>|β|,

    when b<μa, or in other words R0<1. In addition,

    1+γ=2μlμa+μl(μab)<2,

    holds when R0<1. Hence, E0 is locally asymptotically stable when R0<1. In other words, when the larval recruitment is less than adult mortality, the total population will go extinct.

    The positive steady state E1=(L,A), where L,A are defined in (4.1) and (4.2), exists (uniquely) when R0>1. The Jacobian matrix at E1 is given by

    JLA(L,A)=[1μl(1c1A)bec1Aμlec2A1μac2bAe(c1+c2)A]=[1μl(1c1A)bec1Aμlec2A1μac2μaA].

    Here, we have

    β=2μlμac2μaA1+γ=2μaμlc2μaA+μaμlln(R0)=β+μaμlln(R0). (4.6)

    Hence, the first condition for steady state E1 to be asymptotically stable is

    β+μaμlln(bμa)>|β|,

    which is true for R0>1. Hence, |β|<1+γ. According to Lemma 4.2, we also need to provide 1+γ<2. Rewriting Eq (4.6):

    1+γ=2μaμl+μaln(R0)(μlc2c1+c2).

    Since we require 1+γ<2, this implies that

    ln(R0)μa(μlc2c1+c2)<μl+μa.

    For any value of μl>c2c1+c2, we have

    ln(R0)<μa+μlμa(μlc2c1+c2).

    In the case μl<c2c1+c2, this implies that

    ln(R0)>0>μa+μlμa(μlc2c1+c2),

    which reduces to R0>1.

    We present first the global asymptotic stability for the extinction steady sate E0. The following theorem asserts that E0 is globally asymptotically stable when R0<1.

    Theorem 5.1. The steady state E0=(0,0) of system (2.1) and (2.2) is globally asymptotically stable if R0<1. That is, limt(L(t),A(t))=(0,0) when R0<1.

    Proof. Let S(t)=L(t)+A(t), then

    S(t+1)=bA(t)ec1A(t)+(1μl)L(t)+μlL(t)ec2A(t)+(1μa)A(t)bA(t)+(1μl)L(t)+μlL(t)+(1μa)A(t)=S(t)+(bμa)A(t).

    If R0<1, i.e., b<μa, we have S(t+1)S(t). Hence, S(t) is monotone decreasing in t and so limtS(t)=S0. Moreover, since the above inequality implies that

    (μab)A(t)S(t)S(t+1),

    we have by non-negativity that

    lim supt(μab)A(t)limt[S(t)S(t+1)]=limtS(t)limtS(t+1)=0.

    This shows that

    limtA(t)=0.

    From (2.2), we see that

    limtL(t)=limt1μlec2A(t)[A(t+1)(1μa)A(t)]=0.

    To study the global stability of the positive steady state E1=(L,A), we will use the following approach employing a discrete time delay equation from [26]. For simplification, assume sl=1μl and sa=1μa. From Eq (2.2), we have

    L(t1)=11slec2A(t1)[A(t)saA(t1))]. (5.1)

    Using Eqs (2.1), (2.2), and (5.1), we convert system (2.1) and (2.2) into a discrete time delay equation for t2:

    A(t+1)=b(1sl)A(t1)ec1A(t1)c2A(t)+slec2(A(t1)A(t))(A(t)saA(t1))+saA(t). (5.2)

    The dynamics of Eq (5.2) are equivalent to that of system (2.1) and (2.2). The positive steady state of (5.2) is

    ˉA=1c1+c2ln(b1sa),b>1sa,

    which is equivalent to that of the original model (2.1) and (2.2). We note that with this simplified notation, R0=b/(1sa). To study the global stability of ˉA, we use the following theorem from [28] (restated in Theorem 1.1 of [26]).

    Theorem 5.2 ([28]). Consider the difference equation

    xt+1=F(xt,xt1,...,xtk),t0, (5.3)

    where FC(Ω,R), ΩRk+1, is increasing in each of its arguments (when restricted to the region Ω). Let ˉx be an equilibrium of (5.3) and I an interval such that ˉxI and Ik+1={(u1,u2,...,uk+1):uiI,i=1,2,...,k+1}Ω. Moreover, assume that for uˉx,

    (uˉx)(F(u,u,...,u)u)<0.

    Then with initial data x0,x1,...,xkI, we have xtI for t0 and

    limtxt=ˉx.

    Based on the parameter estimation of the LA model, c2 is small (see Table 1). For this reason (and in the spirit of previous analyses [26]), we assume that c2=0 in the rest of this section. For convenience, let xt=A(t+1) for t1. This yields

    xt+1=b(1sl)xt1ec1xt1slsaxt1+(sl+sa)xt, (5.4)

    where,

    x1=A(0),x0=A(1)=(1sl)L(0)+saA(0),

    and the positive steady state of (5.4) is ˉx=1c1ln(b1sa),b>1sa. From Lemma 4.1, we have positive invariance of the solution of system (2.1) and (2.2) in region Ω, which also implies that the solution of Eqs (5.2) and (5.4) is positively invariant to the region when xtbe(1sa)c1. We now establish the global stability of the positive steady state of Eq (5.4) using Theorem 5.2.

    Theorem 5.3. Assume R0>1, so that the positive steady state exists. Assume also that c2=0. If beμlδ, for

    δ=ln(ebμlbμl+e(1μl)(1μa)),

    and 0<b<e, then the positive steady state ˉx of Eq (5.4) with initial data {x1,x0}I, is globally asymptotically stable in the domain Ω=I2, where I=(0,δc1).

    Proof. We consider

    F(xt,xt1)=b(1sl)xt1ec1xt1slsaxt1+(sl+sa)xt, (5.5)

    and

    g(u)=F(u,u)u=b(1sl)uec1u+slsau+sauu.

    According to Lemma 4.1, xt is positive and has an upper bound. Clearly, g(0)=0 and g(+)<0. Hence, for u>0 and, uˉx, we have

    (uˉx)(F(u,u)u)<0.

    Differentiating (5.5),

    Fxt=sa+sl, (5.6)
    Fxt1=b(1sl)(1c1xt1)ec1xt1slsa. (5.7)

    Clearly, F/xt0 as desired. Equation (5.7) may be rewritten as

    Fxt1=b(1sl)ec1xt1c1b(1sl)xt1ec1xt1slsa.

    Recall that xec1x attains a maximum at x=1/c1. Thus,

    Fxt1b(1sl)ec1xt1b(1sl)eslsa.

    Let

    be(1sl)δ=ln(eb(1sl)b(1sl)+eslsa), (5.8)

    then xt1 is positive invariant to region Ω=I2, where I=(0,δc1). By definition, ˉxΩ. When sl,sa(0,1) and b>e, condition (5.8) does not hold. Therefore, for 0<b<e and (5.8), xt1 is restricted to region Ω=I2, ˉxI, and,

    Fxt0andFxt10,when, xt,xt1Ω.

    Therefore, F is increasing in each of its arguments, and we can apply Theorem 5.2. With initial data, x1=A(0)>0 and x0=A(1)=(1sl)L(0)+saA(0)>0, we have limt+xt=ˉx when xt is restricted to Ω. It is worth noting that the conditions presented here are only sufficient for global stability [26].

    Numerically, we can see that in Figure 5, using Parameter Set 1 from Table 1 (c1=0.054,c2=3×104, μl=0.3, and μa=0.2), the solution goes to positive equilibrium state over time for R0>1 and the adult population goes to an extinction state for R0<1. We point out that upper bound of condition 2)(a) of Theorem 4.3, which states that E is locally asymptotically stable in the case μl>c2c1+c2 when

    1<R0<exp(μl+μaμa(μlc2c1+c2))

    can grow arbitrarily large very quickly due to the parameters lying in (0,1). Hence, numerical simulations may appear to indicate global stability at first glance when that is not the case.

    Figure 5.  Solution of (5.2) for the parameter set with c1=0.054,c2=3×104, μl=0.3, μa=0.2, and various b values: b=0.1,0.2,3,10,35,52. Red dotted lines represent the equilibrium points for each value of b. For this set of parameters, b=0.1,0.2 correspond to R0<1, while the other values of b correspond to R0>1.

    We now explore non-steady-state asymptotic behaviors. Varying the adult mortality parameter μa results in potentially chaotic dynamics for approximately μa>0.5, as shown by the chaotic cloud in the left panel and positive Lyapunov exponent in the right panel of Figure 6. This result, obtained using the control group data, is similar to that of the LPA model [3], while the LPAA model in [10] does not exhibit this behavior. The LA model, like the LPA and LPAA models, does not show a period-doubling route to chaos under this parameter regime. This reinforces the idea that model structure informs the complexity exhibited, and may lead to differing results when extrapolating outside of experimental parameters.

    Figure 6.  Bifurcation and Lyapunov exponent diagrams for the LA model fitting with control group data using the parameters b=39.4293, c1=0.0317, c2=4×104, μl=0.4348, and μa(0,1).

    Tribolium populations are often characterized by limit cycles, and some analytical results have been presented [17,29,30]. The LPA model and its variants include transcendental cannibalism terms which complicate analysis. Standard methods for showing the existence and/or stability of a two-cycle involve studying the composition of the model. If p and q with pq are points of a two-cycle for some mapping xt+1=f(xt), then p=f(q) and q=f(p). Then, one could instead look for the fixed points of the composition of the map, i.e., p=f(f(p)). However, the cannibalism terms of the model are transcendental, which limits the use of this technique. Instead, we study the existence of a two-cycle using a numerical approach on the discrete delay form of the LA model (5.2). Let p>0 and q>0 be points of a two-cycle, and pq. For t>1, we have

    p=A(t+1)=A(t1),q=A(t).

    Thus

    p=b(1sl)pec1pc2q+slec2(pq)(qsap)+saq, (6.1)

    and

    q=b(1sl)qec1qc2p+slec2(qp)(psaq)+sap. (6.2)

    A two-cycle exists when (6.1) and (6.2) intersect and pq, although stability is unclear. Some examples are shown in Figure 7. Figure 7(a) shows the curves (6.1) and (6.2) intersecting but with p=q, and hence the simulations do not show a two-cycle. Figure 7(b), (c) show the curves (6.1) and (6.2) intersecting at values where pq, and these intersections give the points of the two cycle. For example, Figure 7(b) has intersections at (p,q)(4,14) and (p,q)(14,4). Due to symmetry, one of these pairs can be considered the points of the two-cycle without loss of generality. While these curves are very dependent on parameter values, no two-cycles are possible when c2=0 due to the only intersection of (6.1) and (6.2) occurring at p=q. Specifically, it appears that a cycle with period nineteen occurs under this parameterization.

    Figure 7.  Implicit plots of (6.1) and (6.2) for b=19, μl=0.54, μa=0.95, and indicated values of c2. The bottom row shows time series simulations for each parameterization. Where all curves intersect, the dashed line indicates a steady-state and is thus not point of the two-cycle.

    In this work, we have built upon and contributed to the study of Tribolium population dynamics by simplifying the LPA model [2] and the LPAA model [10] through the grouping larvae and pupae into one stage (larvae, L), and callow and mature adults into a second stage (adults, A). This decision was driven by new, longitudinal time series for T. confusum in varying nutrient environments in which the populations experienced rapid declines. Unlike the LPA and LPAA models, individuals in the LA model do not undergo continuous, stepwise progression through life stages. Instead, the transition is governed by the unique (1μl) term, which captures the probability of larvae remaining in the same stage rather than advancing to another.

    Several fundamental mathematical properties of the model were established, such as non-negativity and boundedness. The extinction steady state is locally and globally asymptotically stable when R0<1, similar to that of the LPA and LPAA models. As well, for these two models, local stability is particularly sensitive to a ratio of cannibalism intensities (in the case of the LPAA model, adult cannibalism on eggs and callow adults). However, the positive steady state of the LA model is sensitive to the ratio between pupal cannibalism and the sum of the cannibalism intensities. The same technique from Kuang and Cushing [26] can be used to assess the global stability of the positive steady state of the LPAA and LA models. The region for global stability has more constraints than previous work, namely requiring larval recruitment b<e and beμlδ, for

    δ=ln(ebμlbμl+e(1μl)(1μa))

    compared to only constraints on R0 in previous studies. This proof also demonstrates only sufficient conditions for the global stability of the model in the special case c2=0. The LPA and LA models exhibited chaotic dynamics for sufficiently high values of adult mortality μa, while the LPAA model only did for a specific range of values. This reinforces the notion that model choice significantly impacts complexity, particularly when characterizing dynamics outside of experimental regimes. In addition, the existence of two-cycles was explored numerically, illustrating that two-cycles do not exist when cannibalism on pupae c2=0, but appear stable for c2>0.

    This study reinforces that different species of Tribolium and the frequency of media changes impact model performance, previously hypothesized in [10]. However, the poor fitting of the LPAA model may be due to several factors such as the increased ratios of nitrogen and phosphorous added to the media or another biotic or abiotic factor which resulted in declining populations in the experimental and control groups. The differences between the LA and LPAA models are also likely due to the fact that the LA model was optimized for a data set spanning 34 weeks (aggregated data) while the LPAA model was optimized for a data set spanning only 18 weeks.

    Interestingly, this model almost appears continuous in some parameter regimes (see Figure 4 for an example) and demonstrates the potential for a continuous-time model to recapitulate the time-series dynamics. Other formulations of the model could be explored (discrete or continuous time), such as exploring the colonizing and spatial dispersion dynamics of flour beetles or using multi-scale population models. These models (e.g., [31,32]) leverage fractal derivatives to understand interactions at various scales, and may give additional insight into larger patch-type models. Forthcoming work will model Tribolium dynamics using a continuous delay-differential equation system, due to tendency of delay equations to exhibit nonlinear dynamics such as periodic oscillations.

    This research is supported in part by NSF eMB program to SB, YK, JN (DMS-2325146) and IHBEM program DMS-2421258 to YK.

    The authors have no conflicts of interest to declare.

    Code and data availability

    Code and data will be provided upon publication.

    Experimental protocols

    Experimental methods are the same as those provided in Appendix A of [10] and thus summarized here. Longer term studies were conducted of Tribolium responses to variations of inorganic compounds in their environment. A stock colony of T. confusum maintained for the past fifteen years at Scottsdale Community College was used in this study. The population experienced a bottleneck during the COVID-19 pandemic in 2021, reducing to 10% of its maximum. Populations were reestablished after the bottleneck and provided the beetles for this study [10].

    The experiment was started with 20 sclerotized beetles in 25 grams of unbleached flour with 5% yeast media in 118 mL food storage containers at 28C. [10]. This study consisted of several experimental groups where sodium nitrate (NaNO3) and dipotassium phosphate (K2PO4) were added in various ratios. Referring to sodium nitrate as N and dipotassium phosphate as P for brevity, the groups were given by

    1) Control (no nitrogen or phosphorus added)

    2) 5% N

    3) 10% N

    4) 5% P

    5) 10% P

    6) Combination of 2.5% N and 2.5% P (totaling 5%)

    7) Combination of 5% N and 5% P (totaling 10%)

    with percentages by mass. Each group had two trials, with the exception of the control group which had four trials. A larval die off was observed at week 16, potentially due to loss of power.

    Animals were counted every two weeks with media changes every eight weeks. Larvae, pupae, and adults were sifted from the flour and larvae larger than 2.75 mm were considered large enough to count. Eggs are small and blend in with the media, and thus cannot be counted. Callow (unsclerotized)s and dead adults were recorded separately, and if found in pieces dead adults were counted using heads. During media changes, live animals were removed from the media (not including eggs, as they cannot be separated from the media), and 12.5 g of fresh media were added. The number of dead adults removed in the media change was also recorded. For the model, larvae and pupae data were added together, and callow and sclerotized adult data were added together.



    [1] M. D. Pointer, M. J. G. Gage, L. G. Spurgin, Tribolium beetles as a model system in evolution and ecology, Heredity, 126 (2021), 869–883. https://doi.org/10.1038/s41437-021-00420-1 doi: 10.1038/s41437-021-00420-1
    [2] R. F. Costantino, J. M. Cushing, B. Dennis, R. A. Desharnais, Experimentally induced transitions in the dynamic behaviour of insect populations, Nature, 375 (1995), 227–230. https://doi.org/10.1038/375227a0 doi: 10.1038/375227a0
    [3] R. F. Costantino, R. A. Desharnais, J. M. Cushing, B. Dennis, Chaotic dynamics in an insect population, Science, 275 (1997), 389–391. https://doi.org/10.1126/science.275.5298.389 doi: 10.1126/science.275.5298.389
    [4] T. Park, Interspecies competition in populations of Tribolium confusum Duval and Tribolium castaneum Herbst, Ecol. Monogr., 18 (1948), 265–307. https://doi.org/10.2307/1948641 doi: 10.2307/1948641
    [5] T. Park, Experimental studies of interspecies competition Ⅱ. temperature, humidity, and competition in two species of Tribolium, Physiol. Zool., 27 (1954), 177–238. https://doi.org/10.1086/physzool.27.3.30152164 doi: 10.1086/physzool.27.3.30152164
    [6] T. Park, Experimental studies of interspecies competition. Ⅲ relation of initial species proportion to competitive outcome in populations of Tribolium, Physiol. Zool., 30 (1957), 22–40. https://doi.org/10.1086/physzool.30.1.30166306 doi: 10.1086/physzool.30.1.30166306
    [7] M. S. Bartlett, On theoretical models for competitive and predatory biological systems, Biometrika, 44 (1957), 27–42. https://doi.org/10.2307/2333238 doi: 10.2307/2333238
    [8] R. M. May, Biological populations with nonoverlapping generations: Stable points, stable cycles, and chaos, Science, 186 (1974), 645–647. https://doi.org/10.1126/science.186.4164.645 doi: 10.1126/science.186.4164.645
    [9] B. Dennis, R. A. Desharnais, J. M. Cushing, S. M. Henson, R. F. Costantino, Estimating chaos and complex dynamics in an insect population, Ecol. Monogr., 71 (2001), 277–303. https://doi.org/10.1890/0012-9615(2001)071[0277:ECACDI]2.0.CO;2 doi: 10.1890/0012-9615(2001)071[0277:ECACDI]2.0.CO;2
    [10] S. J. Brozak, S. Peralta, T. Phan, J. D. Nagy, Y. Kuang, Dynamics of an LPAA model for Tribolium growth: Insights into population chaos, SIAM J. Appl. Math., 84 (2024), 2300–2320. https://doi.org/10.1137/24M1633881 doi: 10.1137/24M1633881
    [11] H. P. Benoît, E. McCauley, J. R. Post, Testing the demographic consequences of cannibalism in Tribolium confusum, Ecology, 79 (1998), 2839–2851. https://doi.org/10.1890/0012-9658(1998)079[2839:TTDCOC]2.0.CO;2 doi: 10.1890/0012-9658(1998)079[2839:TTDCOC]2.0.CO;2
    [12] T. Park, D. B. Mertz, M. Nathanson, The cannibalism of pupae by adult flour Beetles, Physiol. Zool., 41 (1968), 228–253. https://doi.org/10.1086/physzool.41.2.30155454 doi: 10.1086/physzool.41.2.30155454
    [13] A. Veprauskas, J. M. Cushing, A juvenile–adult population model: Climate change, cannibalism, reproductive synchrony, and strong Allee effects, J. Biol. Dyn., 11 (2017), 1–24. https://doi.org/10.1080/17513758.2015.1131853 doi: 10.1080/17513758.2015.1131853
    [14] S. Via, Cannibalism facilitates the use of a novel environment in the flour beetle, Tribolium castaneum, Heredity, 82 (1999), 267–275. https://doi.org/10.1038/sj.hdy.6884820 doi: 10.1038/sj.hdy.6884820
    [15] J. Cushing, Systems of difference equations and structured population dynamics, in Proceedings of the First International Conference on Difference Equations, (1994), 123–132.
    [16] S. M. Henson, R. F. Costantino, R. A. Desharnais, J. M. Cushing, B. Dennis, Basins of attraction: Population dynamics with two stable 4-cycles, Oikos, 98 (2002), 17–24. https://doi.org/10.1034/j.1600-0706.2002.980102.x doi: 10.1034/j.1600-0706.2002.980102.x
    [17] S. M. Henson, J. M. Cushing, The effect of periodic habitat fluctuations on a nonlinear insect population model, J. Math. Biol., 36 (1997), 201–226, https://doi.org/10.1007/s002850050098 doi: 10.1007/s002850050098
    [18] S. M. Henson, J. M. Cushing, R. F. Costantino, B. Dennis, R. A. Desharnais, Phase switching in population cycles, Proc. R. Soc. London, Ser. B: Biol. Sci., 265 (1998), 2229–2234. https://doi.org/10.1098/rspb.1998.0564 doi: 10.1098/rspb.1998.0564
    [19] S. M. Henson, R. F. Costantino, J. M. Cushing, B. Dennis, R. A. Desharnais, Multiple attractors, saddles, and population dynamics in periodic habitats, Bull. Math. Biol., 61 (1999), 1121–1149. https://doi.org/10.1006/bulm.1999.0136 doi: 10.1006/bulm.1999.0136
    [20] S. M. Henson, Multiple attractors and resonance in periodically forced population models, Physica D: Nonlinear Phenom., 140 (2000), 33–49. https://doi.org/10.1016/S0167-2789(99)00231-6 doi: 10.1016/S0167-2789(99)00231-6
    [21] S. M. Henson, R. F. Costantino, J. M. Cushing, R. A. Desharnais, B. Dennis, A. A. King, Lattice effects observed in chaotic dynamics of experimental populations, Science, 294 (2001), 602–605. https://doi.org/10.1126/science.1063358 doi: 10.1126/science.1063358
    [22] S. M. Henson, J. R. Reilly, S. L. Robertson, M. C. Schu, E. W. Rozier, J. M. Cushing, Predicting irregularities in population cycles, SIAM J. Appl. Dyn. Syst., 2 (2003), 238–253. https://doi.org/10.1137/S1111111102411262 doi: 10.1137/S1111111102411262
    [23] B. Dennis, R. A. Desharnais, J. M. Cushing, R. F. Costantino, Nonlinear demographic dynamics: Mathematical models, statistical methods, and biological experiments, Ecol. Monogr., 65 (1995), 261–282. https://doi.org/10.2307/2937060 doi: 10.2307/2937060
    [24] B. Dennis, R. Desharnais, J. M. Cushing, R. F. Costantino, Transitions in population dynamics: Equilibria to periodic cycles to aperiodic cycles, J. Anim. Ecol., 66 (1997), 704–729. https://doi.org/10.2307/5923 doi: 10.2307/5923
    [25] F. K. Ho, P. S. Dawson, Egg cannibalism by Tribolium larvae, Ecology, 47 (1966), 318–322. https://doi.org/10.2307/1933784 doi: 10.2307/1933784
    [26] Y. Kuang, J. M. Cushing, Global stability in a nonlinear difference-delay equation model of flour beetle population growth, J. Differ. Equ. Appl., 2 (1996), 31–37. https://doi.org/10.1080/10236199608808040 doi: 10.1080/10236199608808040
    [27] L. E. Keshet, Mathematical Models in Biology, SIAM publications, Philadelphia, PA, 2005.
    [28] M. Hautus, T. S. Bolis, A. Emerson, Solution to problem e2721, Am. Math. Mon., 86 (1979), 865–866.
    [29] J. M. Cushing, S. M. Henson, R. A. Desharnais, B. Dennis, R. F. Costantino, A. King, A chaotic attractor in ecology: Theory and experimental data, Chaos Solitons Fractals, 12 (2001), 219–234. https://doi.org/10.1016/S0960-0779(00)00109-0 doi: 10.1016/S0960-0779(00)00109-0
    [30] J. Cushing, Cycle chains and the LPA model, J. Differ. Equ. Appl., 9 (2003), 655–670. https://doi.org/10.1080/1023619021000042216 doi: 10.1080/1023619021000042216
    [31] Y. Zhang, N. Anjum, D. Tian, A. A. Alsolami, Fast and accurate population forecasting with two-scale fractal population dynamics and its application to population economics, Fractals, 32 (2024), 2450082. https://doi.org/10.1142/S0218348X24500828 doi: 10.1142/S0218348X24500828
    [32] N. Anjum, C. H. He, J. H. He, Two-scale fractal theory for the population dynamics, Fractals, 29 (2021), 2150182. https://doi.org/10.1142/S0218348X21501826 doi: 10.1142/S0218348X21501826
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(184) PDF downloads(34) Cited by(0)

Figures and Tables

Figures(7)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog