Is the Primary Helper Always a Key Group for the Sustainability of Cooperative Birds? A Mathematical Study on Cooperative Breeding Blue-Tailed Bee-Eater.

Cooperation is a fundamental requirement for the sustainability of group-living organisms. Despite the substantial research work in cooperative breeding birds, the dependence of the populations’ sustainability young or adult helpers in migratory populations is unidentiﬁed. The mathematical models for predict-ing the birds’ cooperative dynamics so far mostly ig-nore the migratory property. The cooperative breeding birds have three groups in their population, viz, immature or primary helpers, mature or secondary helpers, and breeders. We ask three questions to study migratory cooperative birds’ sustainability through mathematical modeling under changing environments. Which pre-pared group is the key to the sustainability of cooperative migratory birds? Does the maturate helper compensate young helpers’ helping? Does the hierarchical structure of the population vary for variable migratory rates? We explore the answers based on the mathematical model’s simulation experiment, a potential alternative to the game theory approach. This study estimates the parameters associated with the proposed model through the ﬁeld survey and obtains the rest from existing literature. Although the study uses blue-tailed bee-eater as the test-bed species, the model is useful for analyzing other avian species’ behavioral property. The model as a tool can determine whether the primary helpers of blue-tailed bee-eater are the key to sustainability. The model can also classify the adults’ help as an addition or compensate to primary helpers’ help. The model can predict any alteration in the cooperative breeding birds’ hierarchy sizes for variable migration rates under changing climate.

Theoretical Ecology manuscript No. (will be inserted by the editor) group is the key to the sustainability of cooperative migratory birds? Does the maturate helper compensate young helpers' helping? Does the hierarchical structure of the population vary for variable migratory rates? We explore the answers based on the mathematical model's simulation experiment, a potential alternative to the game theory approach. This study estimates the parameters associated with the proposed model through the field survey and obtains the rest from existing literature. Although the study uses blue-tailed bee-eater as the test-bed species, the model is useful for analyzing other avian species' behavioral property. The model as a tool can determine whether the primary helpers of blue-tailed bee-eater are the key to sustainability. The model can also classify the adults' help as an addition or compensate to primary helpers' help. The model can predict any alteration in the cooperative breeding birds' hierarchy sizes for variable migration rates under changing climate.

Introduction
Cooperation between individuals is the yardstick for the population sustainability of group-living animals ranging from insects to birds in the evolutionary timeline (1; 2; 3; 4; 5; 6). The substantial research work on birds' cooperation focuses majorly on cooperative breeding, where non-breeding subordinate helps breeding pairs (7; 8; 9). Ornithologists and behavioral ecologists are interested in finding an association between the bird population sustainability and the helping rate of the subordinate cooperators (10; 11; 12). Identifying the most necessary social-group between breeders and different helpers is another exciting research arena to be explored for behavioral ecologists. Surprisingly, the existing mathematical framework poorly demonstrates the cooperative behavior associated with migration, breeding, and population dynamics. Mathematical models on cooperative breeding determine the behavior's stability and frequency as the opposing force to the Darwinian selection from an evolutionary perspective (12; 7; 13). The cooperative breeding behavior is an essential driver of a group-living organism to eusociality over the evolutionary period (14; 1). So far, the mathematical models enlighten the variation in bird populations' fitness components with a single group of helpers. In general, the research on the cooperative behavior of birds has two domains. The domain of genetic research identifies the genes responsible for inducing cooperation and the helper-breeder relatedness (15; 16). The identification of the cooperative genes is out of the research concern of this article. According to the geneticists' other findings, an animal tends to cooperate with its genetically related kin (17). However, in many birds, non-related adult helpers also cooperate with breeders besides the sexually immature, genetically-related helpers (18; 19; 20). On the other hand, Game theorists focus on the evolutionary stability and repayment of the cooperation against Darwinian selection in avian society. Mathematical ethologists quantified adult cooperators' helping rate for specific environmental conditions (21; 22; 9). Overall the theoretical studies ensure that the adult subordinate birds cooperate as a result of environmental constraints. Their sacrifice in reproductive fitness is repaid almost in the same amount for survival, defense, future breeding opportunities, and social status up-gradation (12; 7; 13). Still, their study does no reveal the population dynamics with variation in the helping rates for both types of helpers under different environments. The output of evolutionary and cooperative game theory proves the helping behavior's evolutionary stability in the social group with multiple strategic situations under natural selection. However, both the sexually mature and immature helpers have the same strategy, i.e., helping breeders to produce and nurture chicks. Suppose two groups of a single population have adopted the same strategy for survival. In that case, the game theory cannot explain further two helpers' role to sustain the population (23). The lacunae in the present theoretical models motivate us to aims at four major questions in cooperative breeding birds.
We use blue-tailed bee-eater (Merops philippinus) as the model species for this study because of their high cooperativeness and long-distance migratory nature. The residence of this bird is far apart from the breeding zone. The model birds' geographic and behavioral distribution allows us to build a model to investigate our objectives in the breeding zone, eliminating residential birds' effect. The other cooperatively breeding birds have different birth and maturation rates. Note that a similar model can understand other birds' behavior with the tuning of the parameter values. The model serves as a general toolbox to study all the cooperative breeding avian societies' ethology when either population size data or the bird's birth rate and maturation rate is known. Commonly, any known cooperative bird's hierarchical classification consists of breeder, primary or immature helpers, and secondary or adult helper. In this present study, we would like to explore three hypotheses: (a) to test whether the sustainability of the model species depends on either primary or secondary helpers; (b) whether the secondary helpers can compensate the primary helpers' cooperation under changing environment if the primary helpers reduce significantly; (c) How does the variation immigration pattern manifest in bird populations' hierarchy sizes? Our theoretical model uses the concept of a nonlinear a delayed differential system equation for breeders, primary helpers, secondary helpers, and prey. We developed a rigorous simulation experiment based on a set of parameter values to obtain this differential equation's solution. These parameter values are not arbitrary, but we obtain most of them from the field survey. We consider most of the remaining parameter values from the existing literature and an online database. The choice of parameter values and its estimation procedure are in the section 3 elaboratively and summarized in table 1.

Model for cooperative breeding birds
We consider three compartments or state variables to represent three social statuses of the cooperatively breeding bird populations for our model. In the model, B(t), H P (t), and H S (t) represent breeders, primary helpers (young chicks), and secondary helpers (adult helpers), respectively. The breeders and the secondary helpers immigrate to the breeding ground at rates i B and i S , respectively. We assume a breeding site can only support B max number of immigrant breeders and H Smax number of adult immigrating helpers due to limited resource. Most long-distance migratory birds, including the blue-tailed bee-eater, migrates asynchronously for the pre-breeding and breeding period (24; 25; 26). Therefore, immigration is not necessarily rapid but sigmoid in most cases. So, we consider a logistic equation to capture the immigrating population. On the other hand, we observed an average rapid emigration of birds from the breeding ground. Hence we denote e B and e S as emigration rates for the B(t) and H ( S), respectively. the equation for immigrant breeders and secondary helpers in the breeding ground is as follows- The breeders (B(t)) produce chicks or primary helpers (H P (t)) at a rate of b in the breeding ground. The chick production rate (b) includes the contributions from other primary and secondary helpers as the sum of cooperation rate per individual helpers (9). If the primary helpers (H P (t)) and the secondary helpers (H S (t)) cooperates at rate h p and h s respectively, the chick production rate b can be expressed as- Here, the b ′ is the chick production rate in absence of helpers. The immature birds leave the breeding ground at a rate of e p and do not come back until they are sexually mature (27). Primary helpers (H P (t)) turn into either breeders (B(t)) or an secondary helpers (H S (t)) upon sexual maturation at respective rate λ b and λ S in the residential zone. Since the maturation takes longer than the migration time, the model has a delay term (τ ) representing the maturation time. Thus the absolute growth rates of breeders, secondary helpers, and primary helpers are- Although the secondary helpers in many birds may breed opportunistically, they do not necessarily form a new breeding pair (7). Already, we captured the rise in chick production rate due to opportunistic breeding in the model parameter h s . Therefore, we do not consider any secondary helpers to breeders turning rate. The birds need one or more prey populations as the food source in both breeding and residential ground. Many cooperative species, especially insectivores, have an extensive range of overlapping prey preferences. So, we consider a nontaxonomic prey population (P (t)). P (t) is an assemble of both migratory and non-migratory insects. Since the insects can come to a breeding site from the nearby sites, we consider a constant introduction rate G P for them. However, the species' death rate (D P ) is constant per individual only due to the predation-pressure by other co-habiting species of the cooperative breeding population. Breeders, primary helpers, and secondary helpers need a different amount of energy as their works in the population are not the same. We consider C B , C p , C S as the breeder's prey consumption rates, primary or young helper, and adult or secondary helper, respectively. The consumed prey will be converted to biomass to sustain the population sizes at different rates too. Let us assume the conversion rate of consumed prey are α B , α p , and α S in order. So, the prey population dynamics will be as follow- Hence, considering the interaction of the birds with prey population the new growth rate of breeders, adult helpers and young helpers or chicks are- respectively. Finally we have the following mathematical model: with the initial biological conditions Remark 1 Primary helpers' helping rate (h p ) involves rates of food sharing, warmth sharing, and parasites' removal resulting in better survival of the chicks. Secondary helpers' helping rate (h s ) is the cumulative measure of the increase in mate finding rate, extra pair fecundity for breeders, rate of food bringing to chicks, and defending chicks from a predator (22). We could not track down the death of birds as it occurs mostly during migration. Hence, the birds' death rates are embedded in the emigration rates.

Model Properties
Our proposed multidimensional model 7 holds non negativity and boundedness (see appendices A and B). The stability analyses of the model find two stable equilibria of the model 7 (see C). One equilibrium is trivial, i.e., no birds exist in the habitat as there is no prey. At the other equilibrium, all populations coexist at different sizes. The theorem, lemma, and proofs for the uniqueness, non-negativity, and equilibrium points are in the appendix 8.
3 Parameter choice using field-data and literature We chose blue-tailed bee-eater as the testbed species. However, the species' population size data is unavailable in avian databases. So, we surveyed seven sites in the breeding grounds of the blue-tailed bee-eaters situated in West Bengal.  figure 1 through "R" software using "mapview" package (28) for visualization of our survey map. At all survey sites, we counted the number of birds coming and leaving per day by point count method. The average of the per hour incoming bird counts are the  Table 1 Collected rate parameters for model sum of immigration rates for both helpers and breeders. We determined the ratio of the secondary helpers and breeders by counting active nest without and with eggs, respectively. The helpers' and breeders' immigration rates were determined by dividing the sum of migration rates by the ratio. We estimated the emigration rates in a similar way from field survey. We observed a breeding pair lay six eggs in an average per nest from our survey. However, all the survey area has adult helpers. Hence the production rate of chicks or the primary helpers from the field survey has already the increment due to helpers. We could not observe the maturation rates of chicks into breeders and helpers in the breeding ground. Also, the rates may vary over different years and different environments (30; 7; 9). So, we could not collect any data on maturation rates from the field survey. The initial values of the rates collected from fields and literature for the model 7 simulation are in the table 1.

Simulation experiment to test the hypotheses
We have simulated the numerical solution for the model through the "DDE23" solver in "MATLAB" software, and parameter values from table 1. Note that the birds' population size differs in various sites. The parameter values may vary for different species too. So our models' numerical solution with table 1 values best describe the blue-tailed bee-eaters' population dynamics in the Gangetic delta plain. Particularly, the maturation time τ varies for different species. It may alter due to the presence of sex-hormone analogs in the environment. As a species-specific property, it determines the time to establish a species' equilibrium in a new habitat. However, we check the models' predictability on other species with different maturation times through simulation. We recorded time to reach the steady-state and size of the steady-state to evaluate the generalization and the species specificity of the model 7.
Under changing environment, the birds' helping rates, maturation rates, and immigration rates may differ. We perform a series of simulation experiments to test our hypotheses for the populations' characteristic under different environmental setups 4.1 Hypothesis 1: Population sustainability depends on the primary and secondary helpers A sustainable population maintains its equilibrium size even under variation in environmental factors. We simulate each populations' equilibrium sizes for a range of helping rates h p and h s . We change one helping rate at a time to test the change in equilibrium size for each helping rate. A paired t-test is performed between the simulated steady-state and the helping rates. The hypothesis is accepted for p − value > 0.05.

Hypothesis 2: Secondary helpers compensate the Primary helpers
The maturation rates into breeder and secondary helpers are both species-specific and environmentally variable. An increased maturation rate from chicks to breeders and secondary helpers can reduce the total help from primary helpers (sibling chicks). We simulate each populations' equilibrium sizes for a range of maturation rates λ b and λ S . We change one maturation rate at a time. A paired t-test is performed between the simulated steady-state and the maturation rates. The hypothesis is accepted for p − value > 0.05.

Hypothesis 3: Populations' hierarchical sizes varies for variable migratory rates
A new habitat may be colonized by a higher number of helpers or a higher number of breeders randomly for a given cooperative bird population. Hence, the immigration rate of breeders and helpers is variable in a new habitat during colonization. We simulate each populations' equilibrium sizes for a range of immigration rates Ii B and i S . A paired t-test is performed between the simulated steady-state and the helping rates. The hypothesis is accepted for p − value > 0.05.

Hierarchy sizes of the population
At our models' coexistence equilibrium, all cooperative breeding bird and prey populations exist together in the habitat. According to our simulated population sizes, the density of the blue-tailed bee-eater's different social hierarchies is different from each other. The primary or the sexually immature helpers are present highest in number, followed by breeders. Adult or secondary helpers are present at the lowest count (figure 2).
Note that the simulated solution has been 100 times downfolded for better visualization. Therefore, the 3.4 number scale of the breeder indicates 340 breeding pairs on average to be expected within 90 days at a random site after the immigration starts (see figure 2). Approximately 750 juvenile helpers can sustain after the first formation of the blue-tailed bee-eater colony. The secondary helper count is around 150, and the prey density is around 7800 per site at the equilibrium. The simulated age structure shows the number of sexually mature birds is less than young birds (B+H S < H P ). Such an age structure is demographically stable in temporal scale (31). The solution also indicates that the breeders and primary helpers achieve a greater bird count value than the equilibrium before reaching the steady-state. However, the cooperative breeding society's secondary or adult helper class decreases rapidly immediately after the first immigration. It then again increases over the period to reach its equilibrium. On the other hand, the birds' prey population shows a rapid increase due to high growth and invasion rate followed by a fall due to the birds' subsequent high prey consumption rate, finally converging into equilibrium. These population dynamics are likely to be observed in emerging colonies of the birds at the breeding zone. The old residential colonies are in equilibrium conditions already.

Species specificity of the model
The variable τ as checked by stability analyses does not affect the equilibrium sizes (see Appendix D). However, the time to achieve a steady-state is indeed affected by variable τ (see figure 3). An extended sexual development period delays the time to reach equilibrium. We find the bird population sizes are large, and the prey population size is low for a species with a short maturation time than blue-tailed bee-eater (see figure 3).

Sustainability of the bird population on different helpers
Our first simulation experiment suggests that the increasing helping rate from sexually immature helpers causes a significant change in the equilibria of all social groups of the bird population and prey size (p-value < 0.05). The higher helping rate from sibling birds raises both the primary helper and secondary helper count. However, it declines the breeder counts and the prey density (See figure 4). Even with the low breeder and prey number, the total bird population keeps increasing till the helping rate is 0.08 per bird (See figure 4). After that prey population drops to zero, and without prey, no birds are found through simulation. Therefore the bird population can not sustain without primary helpers.  Fig. 4 The change in the stable equilibrium sizes of the model (7)  On the contrary, the steady-state of all populations in our model changes insignificantly (pvalue > 0.05) for different helping rates from secondary and adult helpers (see figure 5). We observed almost no changes in the steady-state sizes up to the helping rate of 0.08 per secondary helper, but the helping rate can be further extended. Therefore, the secondary helpers are not necessary for the sustainability of the bird population.

Compensation of Primary helpers' help by Secondary helper
We find a higher value of chicks to breeder maturation rate (λ b ) raises the equilibrium sizes of a breeder, primary helper significantly (p − value < 0.05) in our second simulation experiment. In the same experiment, the prey abundance drops significantly (p−value < 0.05) at the equilibrium. The change in secondary helper count for increasing λ b is insignificant (p − value > 0.05)in our second simulation experiment (see figure 6). So, the secondary helpers do not add a considerable amount of compensation to primary helpers if the breeder population is large. Also, birds will tend to increase the breeder population size given a larger prey population in the environment. Again, for increasing chicks to secondary helpers maturation rate (λ S ), the change in the breeder and prey count's equilibrium size is insignificant in the simulation experiment. However, the primary helper count drops as the secondary helper count increases at the equilibrium in our experiment (see figure 7). Therefore, the secondary helper may compensate primary helpers only if both breeders and primary are limited in the habitat.   Fig. 6 The stable equilibrium sizes of the model (7) populations for variable chicks to breeder maturation rate (λ b ).
(a) Breeder count increases with increasing maturation rate. (b) The chick's count rises with a higher maturation rate.
(c) Secondary helper number increases to a small degree with a higher maturation rate. (d) Prey population size decrease with increased maturation rate from primary helpers to breeder (7) populations.

Effect of migration on helper breeder ratio
Increased breeder immigration decreases the prey population significantly (p−value < 0.05) at the equilibrium in our third simulation experiment. It significantly increases all bird populations at the equilibrium in the experiment (see figure 8). Hence, a change in breeders' migratory pattern can alter the equilibrium sizes of bird populations' hierarchies.
A nine-fold increment of secondary helpers' immigration rate can increase secondary helpers' equilibrium size insignificantly (p−value > 0.05). We observe an insignificant downregulation of The breeder's equilibrium

Discussion
Our models' trivial equilibrium suggested an absence of birds without prey in the habitat. A birds' natural habitat usually have prey. Continuous deforestation and pesticide application cause a high death rate of insects and other prey species. As a consequence, prey species go extinct from some habitats. The cooperative breeding population collapses when the prey population goes extinct in the breeding or residential habitat. One can observe the trivial equilibrium in a habitat before prey comes or after the prey's extinction. A bird will be eliminated from the habitat fast if its prey preference is small. For the equilibrium of coexistence, the model 7 can simulate the bird populations' social hierarchy and portrays the interactions between them. Our simulation is based on blue-tailed bee-eater, but the model is adjustable to all species by putting species-specific maturation time (τ ). The steady-state sizes of the cooperative breeding birds sooner in a habitat with an environment full of sex hormone analogs (e.g., Bisphenol-A for Oestrogen (32; 33)) than other ones. Our proposed model can compare and assess population sustainability's dependency on different helping rates under different environmental conditions. Previous studies on the cooperation rate under different environmental conditions proved that the cooperation increases under restricted resources and greater sex bias (8; 9). Our theoretical study confirms that the primary helpers play an essential to sustain a cooperative breeding bird population. The secondary or adult helpers are insignificant contributors to the sustainability of the avian cooperative breeding system. The increased equilibrium sizes of primary helpers require high consumption of prey to ensure the juveniles' survivability. The secondary helpers can survive on a lower diet. The adult helpers share their prey with juveniles and breeders. In opposition, the breeders require several preys to reproduce more (27). So the number of breeders decreases following the increased primary and secondary helper count in the habitats. Despite a secondary helpers' constant low cooperation, the primary helper can increase their count by helping their siblings and increasing adult helpers. The secondary helpers' opportunistic breeding may be less frequent as their contribution is not essential for the populations' sustainability. The presence of prey species with essential fatty acids in some habitats increases maturation rate into breeder (34). On the other hand, constraints in nesting resources and mates' unavailability force juveniles to mature into non-breeding adult helpers more than into breeders in a sex-biased population. Our model exhibits the same behavior upon simulation. Upon sexual maturation, a high breeder count drops the prey abundance at equilibrium in our model simulation. Our model does not only justify the observations in the existing literature. It also shows that a variable maturation rate into adult helper does not change the prey abundance. So, irrespective of prey abundance, the increased maturation rate can be caused due to species specificity or other demographic effects. For example, catastrophic massdeath occurs under low spatial resources, breeding conditions, and maturation-inducing foods in a small population. Kokko and Johnstone (2001) already stated that a changed ratio of breeders and helpers might occur due to environmental changes. Their study suggested the biased breeder helper ratio increases birds' cooperation by increasing either the helper number or the helping rate.
Their study was specific to a closed population. Hence, the literature so far neglects the bottleneck effect on helper breeder ratio for the migratory population in a given environment. Our simulated bottleneck effect confirms that the helper breeder ratio differs only through a high breeder immigration rate, under two environments for a new habitat. As the population's age structure is sustainable for both in favor of immigration of helpers and breeders, both the helper and breeder may colonize a new habitat at a higher immigration rate. Although the increased breeder immigration increases the total bird population sizes, it may collapse due to a detrimental effect on its prey abundance. Hence a cooperative breeding population must have high breeder immigration. Regardless of the additive or compensatory help from each secondary helper, a high secondary helper immigration rate makes any cooperative breeding population stable in an ecological and evolutionary context. Since different maturation time delay (τ ) does not change populations' overall equilibrium sizes, this model can be generalized for all cooperatively breeding birds. Therefore, the model can be used for classifying the adult helpers' helping rate from any avian population as an additive or compensatory to the juvenile helpers helping rate in a given environment.

Conclusion
The study confirms the cooperative breeding blue-tailed bee-eater population sustain significantly on the primary helpers only. The environmental changes force the secondary helper to increase populations' fitness and their natural fitness. The populations' hierarchy sizes majorly vary for migratory pattern changes. This study has fundamental implications in understanding sustainability through the cooperative behavior of similar species in the avian community.

Declaration
Funding The Council of Scientific and Industrial Research, Human Resource Development, Government of India financially supported the study.

Conflicts of interest
Authors have no conflict of interests.

Availability of data and material
The collected parameter values from field and literature are available in the table 1.

Code availability
The solvers and packages of this study are available online. The authors have customized the codes of models' simulation experiments.

A Non-negativity of solutions
To use the model for a real avian population, the population must survive. The survival of a population in mathematical modelling is possible under the non-negative solution of the model. The results established in (35) for multidimensional models suggests that the solution of our system of equations (7) exists in R 4 + condition. Following the condition, the nonnegativity of the solution can be interpreted and proved using the lemma in (36). In this respect, we have the following theorem.
x 4 (t) = P (t) and f i are right sides of system (7), for example, Now, using the Lemma 1 of (36), we conclude that any solution X(t) of our system equation (7) with X(θ) ∈ C (where, C = ([−τ, 0], R 4 ) is the Banach space of continuous functions with sup-norm), say X(t) = X(t, X(λ)), is such that X(θ) ∈ R 4 + for all t ≥ 0. Therefore, the set R 4 + is an invariant region for system (7).
The existence of non negative solution of our proposed model can now be used for the ethological study if we know the boundedness of our models' solution.

B Boundedness of solution
The boundedness of the solution is useful to depict the model can be described as the region of attraction for the model (7), and provided using the following set: ,ē = min{e P , e B , e S }, m = max{h p , h s }. Each solution that starts in R are ultimately bounded in R. So, in R we can find the equilibrium points of all the populations in this model. The equilibrium points play an important role to under stand the change of sustainability under varying ethological parameters in the model.

C Existence of the equilibrium point
We have checked that the system (7) has one trivial and one equilibrium namely the coexisting equilibrium E * (B * , H * p , H * S , P * ), where (B * , H * p , H * S , P * ) T is a solution of the following system: Since the proposed equation is four dimentional system with multiple interactive mechanism, the steady state values of system (7) is complicated to be determine analytically. A complicated expression of steady states is less interpretable and less efficient as a tool. So, we have determined the steady state numerically while performing the numerical stability analysis on "MATLAB" software.

D Stability of the Equilibrium point
The steady states as checked for our model is either a trivial or a coexistence of all populations. As in trivial points no population exist, the useful steady state from our model is the coexisting one. Still, the model will be only useful if the coexistence equilibrium is stable. For the stability of the equilibrium, we need the Jacobian matrix at the interior equilibrium point at any steady-state. The Jacobian matrix evaluated at the equilibrium point E(B, H P , H S , P ) gives us the following characteristic equation: where I 4 denotes the 4 × 4 identity matrix and M, D, E are defined as follows: where, The characteristic equation can be written as- Since maturation time of a bird cannot be negative, the τ in our model is positive. However,the τ may be too small or approximately 0 in extreme cases of some birds with early maturation due to anthropogenic effects.
D.0.1 Case I: for τ = 0 For τ = 0, the characteristic equation, for τ = 0 at E * , takes the form as Using Routh-Hurwitz (37), we have the following results for the stability of E * for τ = 0.

Proposition 1
The coexistence equilibrium point E * is stable if the following conditions are satisfied: where, D.0.2 Case II: for τ > 0 For τ > 0, (10) will have infinitely many roots. To determine the nature of the stability, the sign of the real parts of the roots of the characteristic equation (10) is required. A necessary condition for a stability changes of E * is that the characteristic equation (10) should have purely imaginary solutions. Let iθ, θ ∈ R, be a root of equation (10). Substituting ρ = iθ in (10) and then separating real and imaginary parts we get then get Let iφ be a root of equation (10). From it we get By squaring and adding the above two equations Substituting φ 2 = ℓ in equation (14), we get the following equation: where, (15) has roots with negative real parts if and only if the Routh-Hurwitz criterion (16) is satisfied. In such case, for these roots we find φ 2 k = ℓ k < 0, so that the roots of (10) are real, iφ k = ∓ √ ℓ k , and not purely imaginary. In summary, we have the following proposition.
Proposition 2 In case of the delayed model (10), the endemic steady state E * is locally asymptotically stable for all τ > 0, if the following conditions are satisfied: Instead, if α 4 < 0 equation (15) possesses at least one positive root m 0 > 0. We find ± i √ m 0 to be a root of (10) corresponding to the delay τ * . By Butler's lemma (38), the endemic equilibrium E * remains stable for τ < τ * , the latter being a critical value of the delay under which the delayed system (7) remains stable. To determine it, from equation 13 we have Thus, in summary we have the following result.
Theorem 2 If α 4 < 0 is satisfied the steady state E * is locally asymptotically stable for τ < τ * and becomes unstable for τ > τ * . Furthermore, the model system will undergo a Hopf-bifurcation at E * when τ = τ * provided that with The first part of the theorem follows from Proposition 3. For the remaining part, we differentiate (10) with respect to τ and get:  (18) and the system undergoes a Hopf bifurcation at τ = τ * . Figure 1 Seven primary data collection sites in breeding zone. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.   The change in the stable equilibrium sizes of the model (7) populations for variable primary helpers' helping rate (hp). (a) Breeder count decreases with increasing helping rate, (b) The increment in helping rate increases raises the primary helpers or chicks count, (c) Secondary helper number increases with higher helping rate, and (d) Prey population size decrease with increased helping ratefrom primary helpers to breeder Figure 5 The change in the stable equilibrium sizes of the model (7)   The stable equilibrium sizes of the model (7) popu-lations for variable chicks to breeder maturation rate (b). (a) Breeder count increases with increasing maturation rate. (b) The chick's count rises with a higher maturation rate. (c) Secondary helper number increases to a small degree with a higher maturation rate.

Figures
(d) Prey population size de-crease with increased maturation rate from primary helpers to breeder (7) populations.

Figure 7
The stable equilibrium sizes of the model (7)    The steady state sizes of the model (7) populations with respect to variable breeder immigration rate. (a) Breeder count falls slightly with rise in breeder immigration rate.(b) Primary helper count increases at a small scale with rise in breeder immigration rate. (c) Secondary helper count in-creases at a small scale with rise in breeder immigration rate. (d) Prey abundance falls slightly with rise in breeder immi-gration rate.