Nutritional regulation influencing colony dynamics and task allocations in social insect colonies

Abstract

In this paper, we use an adaptive modeling framework to model and study how nutritional status (measured by the protein to carbohydrate ratio) may regulate population dynamics and foraging task allocation of social insect colonies. Mathematical analysis of our model shows that both investment to brood rearing and brood nutrition are important for colony survival and dynamics. When division of labour and/or nutrition are in an intermediate value range, the model undergoes a backward bifurcation and creates multiple attractors due to bistability. This bistability implies that there is a threshold population size required for colony survival. When the investment in brood is large enough or nutritional requirements are less strict, the colony tends to survive, otherwise the colony faces collapse. Our model suggests that the needs of colony survival are shaped by the brood survival probability, which requires good nutritional status. As a consequence, better nutritional status can lead to a better survival rate of larvae and thus a larger worker population.

Keywords:

• Social insects
• foraging activities
• nutritional regulation
• backward bifurcation
• bistability dynamics
• adaptive modeling

1. Introduction

In social insect colonies such as ants, bees and wasps, all members of the colony work collectively to ensure colony survival. Colonies act as a single common organism capable of making decisions and forming complex behavioural connections between its members [17,18]. They exhibit a decentralized system with a sophisticated division of labour resulting from interactions among members of the colony and the environment [1,3,18]. In addition to the reproductive division of labour between the queen and the workers, workers also have a division of labour between foragers which leave the nest to search for food and non-foragers that carry out tasks within the nest [9]. In social insect societies, foraging responsibilities are assigned to a subset of adult colony members [8]. Internal and external factors happening at both the individual and colony level shape the foragers' decision to bring back a certain type of food [8].

Currently, there are few studies that have focused on the outcome of nutrient regulation in social insects at the colony level [8] (but see [11,12,32]). Many of these studies lack focus on the overall outcome of colony population dynamics, and how nutrient regulation among foragers affects the number of reared brood, mortality of adult workers, and in general colony survival. In this study, we focus on the mechanisms that regulate foraging behaviour of eusocial workers and the outcomes of these mechanisms on colony performance, including but not limited to the number of brood raised and worker mortality. The collection of food resources by an individual forager is based not only on the colony's current nutritional status but also on the worker's physical caste, age, and prior experience [26,33]. The nutritional needs of the colony are shaped by the differing needs of larvae and workers in the colony [8,10,26]. For instance, the growth of larvae relies heavily on protein, while worker ants require primarily carbohydrates as a source of energy [4,5,10,13,24,34]. Many studies have shown that the ratio of protein to carbohydrates in the diet of a range of insect species is crucial for performance [8,10,21,28,30], though, in general, carbohydrates are often more attractive to foragers than protein [11,26]. However, the protein required for growth may be in greater demand when a queen is laying eggs [26].

In order for social insect foragers to compensate for potential nutrient restrictions in the food available to the colony [6,7,26,29], foragers adjust their collection in favour of food sources containing limiting nutrients [7,11,26]. This guarantees that the colony meets its longer term objectives and thus promotes colony growth and reproduction [26]. According to Dussutour et al. [11], within a colony, workers recruit nestmates for food collection at different rates depending upon food type [5,27], food concentration, and hunger level [11,23]. At the individual level, when workers are starved, recruitment will be stronger to carbohydrate-rich food sources than to sources high in protein [11]. At a collective level, deployment of foragers to carbohydrate-rich or highly proteinaceous material increases in the presence of larvae, resulting in an increase in the collection of carbohydrates and protein [2,11,27].

There are several empirical studies that have studied how a colony is affected by the availability of required nutrients for colony growth and reproduction, and how workers regulate collection of these nutrients to meet individual and collective demands [8,11–13,22,26]. However, currently there are no mathematical models to our knowledge that have attempted to study these mechanisms dynamically. The main goal of this paper is to propose and study an adaptive modeling framework to further understand how nutritional status may regulate population dynamics and foraging task allocation of social insect colonies. The proposed model contains three compartments that allow us to analyse and measure the impacts of nutritional status that can benefit colony growth and survival. Our model assumes that (1) nutritional status is measured by the protein to carbohydrate ratio, which reflects the ratio of workers foraging for protein to workers foraging for carbohydrates; (2) brood are able to survive if the protein to carbohydrate ratio falls into a certain range; and (3) the colony recruits workers to forage for protein or carbohydrate in order to maximize the brood survival rate. In addition, our proposed model includes division of labour implicitly. Also, by considering the basic mechanisms affecting colony growth such as cooperative effort for reproductive division of labour, successful brood maturation/survival, and recruitment of workers to collect different nutrients based on specific colony nutritional demands, our model could help us understand how other life history factors affect the performance (number of brood raised and mortality of workers) of the colony.

The rest of this article is organized as follows: In Section 2, we describe the detailed derivation of our proposed model. In Section 3, we provide the mathematical analysis of our model including lemmas, propositions, and theorems, the proofs of which can be found in Section 6. In Section 4, we provide numerical simulations illustrating the equilibrium dynamics of the model to further obtain biological insights of some life history parameters of the colony. Lastly, the conclusion of this paper is found in Section 5.

2. Derivation of the mathematical model

Let $L(t)$ represent the brood population; $A_p(t)$ be the portion of foragers collecting proteinaceous material, called the protein forager; $A_c(t)$ be the portion of foragers collecting carbohydrates, called the carbohydrate forager. The total forager population is denoted as $A(t)=A_p(t)+A_c(t)$. The following ecological assumptions determine the population of $L,A_c$ and $A_p$:

1. Brood population $L(t)$: The brood population L increases with the average egg-laying rate of the queen(s) given by γ, which is discounted by two factors:

   1. (a) The survival rate function of eggs is determined by the cooperative efforts of workers A in the colony. We adopt the modeling approach from Kang et al. [15,25], where the cooperative efforts that lead to the eggs' survival is measured by a Holling type-III function $\frac{a{A}^2}{b+a{A}^2}$, where b is a half-saturation constant and a is the portion of the division of labour invested towards the successful development of the larvae.

   2. (b) The survival rate of larvae to workers is determined by the available nutrients in the colony which is reflected through the protein to carbohydrate ratio of worker collectors $S_L\left(\frac{A_p}{A_c}\right)$. Examples of $S_L$ could be $S_L\left(\frac{A_p}{A_c}\right)=-{\alpha }_1|\frac{A_p}{A_c}-{\theta }_m|+{\alpha }_2({\theta }_c-{\theta }_m)$ with ${\alpha }_i\in (0,1)i=1,2$ as a scaling factor of nutrient collection, ${\theta }_m$ representing the optimal nutrient ratio, and ${\theta }_c$ representing the maximal nutrient ratio that brood can survive (see Figure 1(a)), or general functions such as the normal biological performance curve (see Figure 1(b)). Notice that $S_L\le 1$ can be negative, thus we define $S_{L_{max}}=max\{0,S_L\}$ such that $S_{L_{max}}\in [0,1]$ is a survival probability.

   The brood population decreases by a maturation rate $\beta L$, which describes the rate at which brood matures into the adult class A. Thus, we have following equation: $L^{\prime }=\gamma \cdot \underset{\text{nutrient effects}}{\underbrace{S_{L_{max}}}}\cdot \underset{\text{adult worker efforts}}{\underbrace{\frac{a{A}^2}{b+a{A}^2}}}-\underset{\text{maturation rate}}{\underbrace{\beta L}}.$ When $\frac{A_p}{A_c}$ is less than ${\theta }_m$, the brood survival rate increases, and decreases when $\frac{A_p}{A_c}$ is greater than ${\theta }_m$. This phenomenon has been supported by the work of [8,12,20,21], in which it is explained that worker survivability decreases as a probable side effect of an excess ingestion of proteins and of carbohydrate limitation. Figure 1(a) shows a general case of $S_L\left(\frac{A_p}{A_c}\right)=-{\alpha }_1|\frac{A_p}{A_c}-{\theta }_m|+{\alpha }_2({\theta }_c-{\theta }_m)$ with different ${\alpha }_1$ and ${\alpha }_2$. The partial derivative of $S_L=S_L\left(\frac{A_p}{A_c}\right)$ reveals constant rates, showing a linear relation between the brood survival rate and nutritional status. In this study, we assume that when the nutrition level hits or exceeds the critical value ${\theta }_c$, i.e. $\frac{A_p}{A_c}\ge {\theta }_c$, the nutrient becomes toxic such that no brood can survive. Lastly, Figure 1(b) shows how the survival rate grows with respect to the collection of nutrients until it reaches ${\theta }_m$.

   In general, we expect the protein to carbohydrate ratio $\frac{A_p}{A_c}$ of the colony to fall in a certain range in order for the colony to survive and grow, say, $\frac{A_p}{A_c}\in [{\theta }_0,{\theta }_m]$, where ${\theta }_0$ is the minimum nutrient ratio necessary for brood survival. This is supported by [8,12]. Thus it is reasonable to assume that $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$ has the following simple form: (1) $S_{L_{max}}\left(\frac{A_p}{A_c}\right)=\{\begin{array}{ll}0& \text{ when }0\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_0,\\ {\alpha }_1({\displaystyle \frac{A_p}{A_c}}-{\theta }_0)& \text{when }{\theta }_0\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_m,\\ -{\alpha }_2({\displaystyle \frac{A_p}{A_c}}-{\theta }_c)& \text{when }{\theta }_m\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_c,\\ 0& \text{when }{\theta }_c\le {\displaystyle \frac{A_p}{A_c}},\end{array}$(1) subject to ${\alpha }_i\in (0,1),i=1,2$ and ${\alpha }_1({\theta }_m-{\theta }_0)+{\alpha }_2({\theta }_m-{\theta }_c)=0\text{ and }0<{\alpha }_1({\theta }_m-{\theta }_0)\le 1.$

   In particular, we have $S_L={\alpha }_1(\frac{A_p}{A_c}-{\theta }_0)\text{ when, }\frac{A_p}{A_c}<{\theta }_m;S_L=-{\alpha }_2(\frac{A_p}{A_c}-{\theta }_c)\text{ when, }\frac{A_p}{A_c}\ge {\theta }_m.$ Then we have (2) $\begin{array}{rl}& \frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}={\alpha }_1;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}=-{\alpha }_1{\left(\frac{A_p}{A_c}\right)}^2<0\text{when }\frac{A_p}{A_c}<{\theta }_m,\\ & \frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}=-{\alpha }_2;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}={\alpha }_2{\left(\frac{A_p}{A_c}\right)}^2>0\text{when }{\theta }_m\le \frac{A_p}{A_c}<{\theta }_c.\end{array}$(2) In the symmetric case, i.e. ${\alpha }_1={\alpha }_2=\alpha$, then we have $S_L\left(\frac{A_p}{A_c}\right)=-\alpha |\frac{A_p}{A_c}-{\theta }_m|+\alpha ({\theta }_c-{\theta }_m)\text{when }\frac{A_p}{A_c}\in [{\theta }_0,{\theta }_c]$ with $\alpha \in (0,1)$, ${\theta }_0=2{\theta }_m-{\theta }_c$ and ${\theta }_c\in [{\theta }_m,2{\theta }_m+\frac{1}{\alpha }]$. Thus, $S_L=\alpha (\frac{A_p}{A_c}+{\theta }_c-2{\theta }_m)$ if $2{\theta }_m-{\theta }_c\le \frac{A_p}{A_c}<{\theta }_m$ and $S_L=\alpha ({\theta }_c-\frac{A_p}{A_c})$ if ${\theta }_m\le \frac{A_p}{A_c}<{\theta }_c$. In addition, we have the following: (3) $\begin{array}{rl}& \frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}=\alpha ;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}=-\alpha {\left(\frac{A_p}{A_c}\right)}^2\text{ ,when, }2{\theta }_m-{\theta }_c\le \frac{A_p}{A_c}<{\theta }_m,\\ & \frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}=-\alpha ;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}=\alpha {\left(\frac{A_p}{A_c}\right)}^2\text{ ,when, }{\theta }_m\le \frac{A_p}{A_c}<{\theta }_c.\end{array}$(3) The special case of the symmetric scenario above is $S_L(0)=0$, i.e. ${\theta }_c=2{\theta }_m$. In this case, we have $S_L\left(\frac{A_p}{A_c}\right)=-\alpha |\frac{A_p}{A_c}-{\theta }_m|+\alpha {\theta }_m$ with ${\theta }_m\in (0,\frac{1}{\alpha })$ in our proposed model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ). In the following section, we will provide mathematical analysis of the general case of $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$ shown in (1(1) $S_{L_{max}}\left(\frac{A_p}{A_c}\right)=\{\begin{array}{ll}0& \text{ when }0\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_0,\\ {\alpha }_1({\displaystyle \frac{A_p}{A_c}}-{\theta }_0)& \text{when }{\theta }_0\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_m,\\ -{\alpha }_2({\displaystyle \frac{A_p}{A_c}}-{\theta }_c)& \text{when }{\theta }_m\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_c,\\ 0& \text{when }{\theta }_c\le {\displaystyle \frac{A_p}{A_c}},\end{array}$(1) ) and the related results can be applied to the symmetric case and its special case directly.

2. The total forager population $A(t)$: The population A increases by the maturation rate of brood and decreases with a density-dependent death rate $d{A}^2$. The density-dependent mortality rate follows the approach of [15], thus we have following equation: $A^{\prime }=\underset{\text{maturation from brood}}{\underbrace{\beta L}}-\underset{\text{average mortality rate}}{\underbrace{d{A}^2}}.$

3. The protein forager population $A_p(t)$: The ratio $\frac{A_p}{A_c}$ measures the nutritional status of the colony. Assume that the brood can survive in a range of nutrient ratio, i.e. $\frac{A_p}{A_c}\in ({\theta }_0,{\theta }_c)$, and any ratio greater than ${\theta }_c$ can be toxic to brood. In addition, there is an optimal nutritional ratio $\frac{A_p}{A_c}$, denoted by ${\theta }_m$, such that brood could have the optimum survival rate at this ratio. More specifically, the brood survival rate increases with respect to the value of $\frac{A_p}{A_c}$ when $\frac{A_p}{A_c}\in ({\theta }_0,{\theta }_m)$, and passing this optimal ratio ${\theta }_m$, the brood survival rate decreases with $\frac{A_p}{A_c}$. The survival rate of brood is zero when $\frac{A_p}{A_c}\le {\theta }_0$ or $\frac{A_p}{A_c}\ge {\theta }_c$. Thus,

   	•	The portion of the successful brood developed into adults which enter into the protein forager pool can be modeled by the term: $\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\},$ where $max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}\in [0,1]$ represents the nutritional requirements of the colony measured by the ratio $\frac{A_p}{A_c}$ and a nutrient collection factor.
   	•	Based on the nutritional requirements of the colony and other related stimuli, a protein forager can become a carbohydrate forager, and vice versa. This task switching rate depends upon different factors such as the nutritional status of the colony, presence of larvae, individual preference, food type, food concentration and hunger level [5,7,11,23,26,27]. In this paper, we assume that the task switching rate of protein foragers to carbohydrate foragers depends on the brood population L, the nutritional requirement of the colony $max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}$, and the available carbohydrate forager $A_c=A-A_p$, thus its switching rate is $max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_{p}L$. Similarly, the switching rate of carbohydrate foragers to protein foragers is termed as $max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_{c}L$. This gives the net task switching rate of the protein forager: $(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L$. For instance, if the ratio $\frac{A_p}{A_c}$ is less than the optimal nutrient ratio ${\theta }_m$ (where the maximum brood survival rate occurs), then we expect $\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}>0$ and $\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}<0$, thus this indicates that carbohydrate foragers will switch tasks to forage for protein, i.e. $max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_{p}L=0$. In a similar fashion, if the ratio ${\theta }_m\le \frac{A_p}{A_c}\le {\theta }_c$, protein foragers will switch to forage for carbohydrates, that is, $max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_{c}L=0$.
   	•	The total forager population $A=A_p+A_c$ decreases with a density-dependent death rate $d{A}^2=dA(A_p+A_c)$, then the protein forager population decreases with the density-dependent mortality rate $dA{A}_p$.

   Considering the factors above, we derive the population dynamics of the protein forager as follows: $\begin{array}{rl}{A}_p^{\prime }& =\underset{\text{portion of matured adults entering {}{A}_p\text{}}}{\underbrace{\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}}}\\ & +\underset{\text{net task switching}}{\underbrace{(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L}}-\underset{\text{mortality rate}}{\underbrace{dA{A}_p}}.\end{array}$

Figure 1. Examples of a survival rate function: (a) A general case of $S_L\left(\frac{A_p}{A_c}\right)=-{\alpha }_1|\frac{A_p}{A_c}-{\theta }_m|+{\alpha }_2({\theta }_c-{\theta }_m)$ with different ${\alpha }_1=0.3$ and ${\alpha }_2=0.15$; (b) The normal biological performance curve $S_L\left(\frac{A_p}{A_c}\right)$.

Based on the ecological assumptions above, the population dynamics of a social insect colony with nutrient regulating foraging activities is described as follows: (4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) The biological meaning of the parameters and the related values are listed in Table 1.

Table 1. Parameter description and interval values of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ).

Parameter	Description	Range	Reference
a	Portion of the division of labour invested on larvae	$(0,0.25)$	[15]
b	Half-saturation constant	$(0.001,10)$	[15]
d	Adult worker death rate	$(0.001,1)$	[15]
$\alpha ,{\alpha }_i,i=1,2$	Scaling factor(s) of nutrient collection	$(0,1)$
β	Maturation rate from larvae to adult worker	$(0,1)$	[7]
γ	Average egg-laying rate of queen	$(0,1)$	[11]
${\theta }_m$	Optimal nutrient ratio	$(0,{\theta }_c)$
${\theta }_c$	Maximal nutrient ratio that brood can survive	$(0,\frac{1}{\alpha }+{\theta }_m)$
${\theta }_0$	Minimal nutrient ratio for brood survival	$(0,{\theta }_m)$

3. Mathematical analysis

The state space of the proposed ecological model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is ${\mathbb{R}}_{+}^3$. All parameters $a,b,d,\alpha ,\beta ,\gamma$, ${\theta }_0,{\theta }_m,{\theta }_c$ are assumed to be strictly positive based on their biological meaning. We focus on the proposed function $S_L\left(\frac{A_p}{A_c}\right)$ shown in Equation (1(1) $S_{L_{max}}\left(\frac{A_p}{A_c}\right)=\{\begin{array}{ll}0& \text{ when }0\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_0,\\ {\alpha }_1({\displaystyle \frac{A_p}{A_c}}-{\theta }_0)& \text{when }{\theta }_0\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_m,\\ -{\alpha }_2({\displaystyle \frac{A_p}{A_c}}-{\theta }_c)& \text{when }{\theta }_m\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_c,\\ 0& \text{when }{\theta }_c\le {\displaystyle \frac{A_p}{A_c}},\end{array}$(1) ) and Figure 1(a). The related mathematical results should be easily adopted to the symmetric case (3(3) $\begin{array}{rl}& \frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}=\alpha ;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}=-\alpha {\left(\frac{A_p}{A_c}\right)}^2\text{ ,when, }2{\theta }_m-{\theta }_c\le \frac{A_p}{A_c}<{\theta }_m,\\ & \frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}=-\alpha ;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}=\alpha {\left(\frac{A_p}{A_c}\right)}^2\text{ ,when, }{\theta }_m\le \frac{A_p}{A_c}<{\theta }_c.\end{array}$(3) ). Under such conditions, we first show that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is biologically well-defined, i.e. it is positively invariant and bounded in ${\mathbb{R}}_{+}^3$ in the following lemma:

Lemma 3.1

Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is positively invariant and bounded in ${\mathbb{R}}_{+}^3=\{(L,A,A_p):L\ge 0,A\ge 0,A_p\ge 0\}$. In particular, if $L(0)>0,A(0)>0$ and $A_p(0)>0,$ then $L(t)>0,A(t)>0$ and $A_p(t)>0$ for all t>0.

The extinction equilibrium $E_0=(0,0,0)$ of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) always exists. The local stability of the trivial equilibrium $E_0$ cannot be analysed directly for our model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ). However, from the first two equations of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ), we have $\begin{array}{rl}(L+A{)}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-d{A}^2\le [\frac{\gamma {\alpha }_1({\theta }_m-{\theta }_0)a}{b+a{A}^2}-d]A^2\\ & \le [\frac{\gamma {\alpha }_1({\theta }_m-{\theta }_0)a}{b}-d]A^2.\end{array}$ Note that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is positively invariant and bounded from Lemma 3.1, thus we can conclude that if $\frac{a{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}{b}<d$, then the inequality above implies that $\underset{t\to \mathrm{\infty }}{lim sup}(L+A)$ converges to a nonnegative constant. In addition, we have $L^{\prime }{|}_{A=0,L>0}=-\beta L<0$, $A^{\prime }{|}_{L=0}=-d{A}^2<0$ and ${A_p}^{\prime }{|}_{L=0,A>0}=-d{A}_{p}A<0$. Therefore, if $\frac{a{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}{b}<d$, then for some initial conditions around $E_0=(0,0,0)$, Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) converges to the extinction equilibrium $E_0$ in ${\mathbb{R}}_{+}^3$. Thus, we have the following proposition:

Proposition 3.1

If $\frac{a{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}{b}<d,$ then for some initial condition around the extinction equilibrium point $E_0=(0,0,0),$ taken in ${\mathbb{R}}_{+}^3,$ the trajectory of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) converges to $E_0$.

Remark 3.1

Note that the inequality $\frac{a{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}{b}<d$ implies that $a<\frac{bd}{\gamma {\alpha }_1({\theta }_m-{\theta }_0)}$. Proposition 3.1 implies that if a is not large enough (i.e. the investment to the brood growth is small), or the death rate of adults is too large, then the brood population and the total forager population approaches the extinction equilibrium point $E_0$. In the symmetric case, we have ${\theta }_0=2{\theta }_m-{\theta }_c$, then the inequality becomes $\frac{a\alpha \gamma ({\theta }_c-{\theta }_m)}{b}<d$ with ${\alpha }_1={\alpha }_2=\alpha$.

Assume that $E^{\ast }=(L^{\ast },A^{\ast },A_p^{\ast })$ is an interior equilibrium of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) with the general case of $S_L$. Then based on the equation of $\frac{\mathrm{d}{A}_p}{\mathrm{d}t}$ shown in (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ), we can conclude that $A_p^{\ast }$ can exist only if $\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}>0$ as it requires $max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}>0$. Biologically, this implies that the colony survival requires the nutritional needs of brood being on the positive gradient of the brood survival probability $S_{L_{max}}$. Thus, we have $\frac{A_p^{\ast }}{A_c^{\ast }}=\frac{A_p^{\ast }}{A^{\ast }-A_p^{\ast }}\in ({\theta }_0,{\theta }_m)$, and therefore $S_L\left(\frac{A_p}{A_c}\right)={\alpha }_1(\frac{A_p}{A_c}-{\theta }_0)$ and $\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}{|{}_{A_c=A_c^{\ast },A_p=A_p^{\ast }}={\alpha }_1;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}|}_{A_c=A_c^{\ast },A_p=A_p^{\ast }}=-{\alpha }_1{\left(\frac{A_p^{\ast }}{A_c^{\ast }}\right)}^2<0.$ To solve for $(L^{\ast },A^{\ast },A_p^{\ast })$, we set $L^{\prime }=A^{\prime }=A_p^{\prime }=0$, which implies the following equations $\begin{array}{rl}{L}^{\prime }& =0⟹{\alpha }_1\gamma [\frac{A_p}{A_c}-{\theta }_0]\frac{a{A}^2}{b+a{A}^2}-\beta L=0,\\ A^{\prime }& =0⟹L=\frac{d}{\beta }{A}^2,\\ A_p^{\prime }& =0⟹{\alpha }_1\beta L+{\alpha }_1{A}_{c}L-d{A}_{p}A=0⟹L=\frac{d{A}_{p}A}{{\alpha }_1\beta +{\alpha }_1{A}_c},\end{array}$ which gives $L=\frac{d}{\beta }{A}^2\text{and}L=\frac{d{A}_{p}A}{{\alpha }_1\beta +{\alpha }_1{A}_c}.$ Therefore, $A^{\ast }$ of an interior equilibrium $(L^{\ast },A^{\ast },A_p^{\ast })$ satisfies the following equation: (5) $ad\beta (1-{\alpha }_1)A^2-a{\alpha }_1^2\gamma A+\beta [(1-{\alpha }_1)(bd+a{\alpha }_1\gamma {\theta }_0)-a{\alpha }_1^2\gamma ]=0.$(5) Recall that ${\alpha }_1\in (0,1)$. Depending on the exact values of $a,b,d,{\alpha }_1,\beta ,\gamma ,{\theta }_0,{\theta }_m,{\theta }_c$, the Equation (5(5) $ad\beta (1-{\alpha }_1)A^2-a{\alpha }_1^2\gamma A+\beta [(1-{\alpha }_1)(bd+a{\alpha }_1\gamma {\theta }_0)-a{\alpha }_1^2\gamma ]=0.$(5) ) can have zero, one, or two positive roots.

Let $A_1=\frac{a{\alpha }_1^2\gamma -\sqrt{\mathrm{\Delta }}}{2a\beta d(1-{\alpha }_1)},A_2=\frac{a{\alpha }_1^2\gamma +\sqrt{\mathrm{\Delta }}}{2a\beta d(1-{\alpha }_1)},$ where $\mathrm{\Delta }=a(a{\alpha }_1^4{\gamma }^2-4d{\beta }^2[(1-{\alpha }_1{)}^2(bd+a{\alpha }_1\gamma {\theta }_0)-a{\alpha }_1^2\gamma (1-{\alpha }_1)])$ be the possible positive roots of Equation (5(5) $ad\beta (1-{\alpha }_1)A^2-a{\alpha }_1^2\gamma A+\beta [(1-{\alpha }_1)(bd+a{\alpha }_1\gamma {\theta }_0)-a{\alpha }_1^2\gamma ]=0.$(5) ). Let us denote ${\hat{a}}^{\ast }$ as follows: (6) $\begin{array}{rl}{\hat{a}}^{\ast }& =\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^4{\gamma }^2+4d{\beta }^2{\alpha }_1\gamma (1-{\alpha }_1)[{\alpha }_1-{\theta }_0(1-{\alpha }_1)]}\\ & =\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}\frac{4d{\beta }^2(1-{\alpha }_1)}{\frac{{\alpha }_1^3\gamma }{{\alpha }_1-(1-{\alpha }_1){\theta }_0}+4d{\beta }^2(1-{\alpha }_1)}<\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}\end{array}$(6) which is an increasing function of ${\theta }_0$.

In the symmetric case, we have ${\theta }_0=2{\theta }_m-{\theta }_c$, then ${\hat{a}}^{\ast }$ shown in (6(6) $\begin{array}{rl}{\hat{a}}^{\ast }& =\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^4{\gamma }^2+4d{\beta }^2{\alpha }_1\gamma (1-{\alpha }_1)[{\alpha }_1-{\theta }_0(1-{\alpha }_1)]}\\ & =\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}\frac{4d{\beta }^2(1-{\alpha }_1)}{\frac{{\alpha }_1^3\gamma }{{\alpha }_1-(1-{\alpha }_1){\theta }_0}+4d{\beta }^2(1-{\alpha }_1)}<\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}\end{array}$(6) ) can be rewritten as (7) $\begin{array}{r}{\tilde{a}}^{\ast }=\frac{4b{d}^2{\beta }^2(1-\alpha {)}^2}{{\alpha }^4{\gamma }^2+4d\alpha {\beta }^2\gamma (1-\alpha )[\alpha -(1-\alpha )(2{\theta }_m-{\theta }_c)]}<\frac{bd(1-\alpha )}{\alpha \gamma [\alpha -(1-\alpha )(2{\theta }_m-{\theta }_c)]}.\end{array}$(7) Also note that $\frac{{\alpha }_1^3\gamma +4d{\alpha }_1{\beta }^2(1-{\alpha }_1)}{4d{\beta }^2(1-{\alpha }_1{)}^2}=\frac{{\alpha }_1^3\gamma }{4d{\beta }^2(1-{\alpha }_1{)}^2}+\frac{{\alpha }_1}{1-{\alpha }_1}>\frac{{\alpha }_1}{1-{\alpha }_1}$. Then the following theorem provide conditions for existence of equilibrium solutions of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ):

Theorem 3.1

Existence of Equilibria

For Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ),

1.	If $0<a<{\hat{a}}^{\ast }$ and ${\theta }_0<\frac{{\alpha }_1^3\gamma +4d{\alpha }_1{\beta }^2(1-{\alpha }_1)}{4d{\beta }^2(1-{\alpha }_1{)}^2}$, then there is only one trivial equilibrium $E_0=(0,0,0)$ and no other positive equilibrium.
2.	If $a={\hat{a}}^{\ast }$ and ${\theta }_0<\frac{{\alpha }_1^3\gamma +4d{\alpha }_1{\beta }^2(1-{\alpha }_1)}{4d{\beta }^2(1-{\alpha }_1{)}^2}$, then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two positive equilibria which collapse into one equilibrium $E_{\ast }$ $E_{\ast }=(L_{\ast },A_{\ast },A{}_{p_{\ast }})=(\frac{d}{\beta }{A}_{\ast }^2,\frac{{\alpha }_1^2\gamma }{2\beta d(1-{\alpha }_1)},\frac{{\alpha }_1\beta A_{\ast }+{\alpha }_1{A}_{\ast }^2}{\beta +{\alpha }_1{A}_{\ast }})$ in addition to $E_0=(0,0,0)$.
3.	If $a>\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}$ and ${\theta }_0<\frac{{\alpha }_1}{1-{\alpha }_1}$, then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has only one positive equilibrium $E_2$ $E_2=(\frac{d}{\beta }{A}_2^2,A_2,\frac{{\alpha }_1\beta A_2+{\alpha }_1{A}_2^2}{\beta +{\alpha }_1{A}_2})$ in addition to $E_0$.
4.	If ${\hat{a}}^{\ast }<a<\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}$ and ${\theta }_0<\frac{{\alpha }_1}{1-{\alpha }_1}$, then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two positive equilibria in the following form in addition to $E_0$: $E_1=(\frac{d}{\beta }{A}_1^2,A_1,\frac{{\alpha }_1\beta A_1+{\alpha }_1{A}_1^2}{\beta +{\alpha }_1{A}_1})\mathrm{and}{E}_2=(\frac{d}{\beta }{A}_2^2,A_2,\frac{{\alpha }_1\beta A_2+{\alpha }_1{A}_2^2}{\beta +{\alpha }_1{A}_2}).$

Remark 3.2

The detailed proof of Theorem 3.1 is shown in the last section. The number of equilibria of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is determined by the positive root(s) of Equation (5(5) $ad\beta (1-{\alpha }_1)A^2-a{\alpha }_1^2\gamma A+\beta [(1-{\alpha }_1)(bd+a{\alpha }_1\gamma {\theta }_0)-a{\alpha }_1^2\gamma ]=0.$(5) ). Theorem 3.1 implies that the value of the division of labour invested on larvae a and the minimal protein to carbohydrate ratio ${\theta }_0$ determine the existence of the interior equilibrium $(L_i,A_i,A_{pi}),i=1,2$.

Our simulations (see Section 4) suggest that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has simple dynamics: no limit cycle and only equilibrium dynamics. At the stable equilibrium, the ratio describing the nutritional level of the colony is (8) $\frac{A_p^{\ast }}{A_c^{\ast }}=\frac{{\alpha }_1(A^{\ast }+\beta )}{\beta (1-{\alpha }_1)}\in ({\theta }_0,{\theta }_c).$(8) Equation (8(8) $\frac{A_p^{\ast }}{A_c^{\ast }}=\frac{{\alpha }_1(A^{\ast }+\beta )}{\beta (1-{\alpha }_1)}\in ({\theta }_0,{\theta }_c).$(8) ) suggests that the larger the total population of workers investing in nutrient collection is, the higher the ratio of protein to carbohydrates will be, i.e. better nutrient status of the colony. Notice that $A^{\ast }$ depends on ${\theta }_0$, so $\frac{A_p^{\ast }}{A_c^{\ast }}$ does as well.

Now we discuss stability of the interior equilibrium for Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ). Let $E^{\ast }=(L^{\ast },A^{\ast },A_p^{\ast })$ be an arbitrary positive interior equilibrium of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ). The Jacobian matrix associated to Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) at equilibrium is: (9) $\begin{array}{rl}J{|}_{E^{\ast }}& =\left(\begin{array}{ccc}-\beta & J_{12}& J_{13}\\ \beta & -2d{A}^{\ast }& 0\\ J_{31}& J_{32}& J_{33}\end{array}\right),\\ J_{12}& =\frac{a{\alpha }_1\gamma A^{\ast }[A_p^{\ast }(b{A}^{\ast }-a{A}^{\ast 3}-2b{A}_p^{\ast })-2b(A^{\ast }-A_p^{\ast }{)}^2{\theta }_0]}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A}^{\ast 2}{)}^2},\\ J_{13}& =\frac{a{\alpha }_1\gamma A^{\ast 3}}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A}^{\ast 2})}>0,\\ J_{31}& ={\alpha }_1(\beta +A^{\ast }-A_p^{\ast })>0,J_{32}={\alpha }_1{L}^{\ast }-d{A}_p^{\ast },J_{33}=-({\alpha }_1{L}^{\ast }+d{A}^{\ast })<0.\end{array}$(9) Then the characteristic equation of $J{|}_{E^{\ast }}$ is (10) $f(\lambda )={\lambda }^3+C_1{\lambda }^2+C_2\lambda +C_3=0,$(10) where (11) $\begin{array}{rl}{C}_1& =\beta +\alpha L^{\ast }+3d{A}^{\ast }>0,\\ C_2& =J_{11}{J}_{33}+J_{11}{J}_{22}+J_{22}{J}_{33}-J_{21}{J}_{12}-J_{31}{J}_{13},\\ C_3& =-\det (J{|}_{E_i^{\ast }})=J_{11}{J}_{22}{J}_{33}+J_{21}{J}_{32}{J}_{13}-J_{21}{J}_{12}{J}_{33}-J_{31}{J}_{22}{J}_{13}.\end{array}$(11) The stability of the steady state $E^{\ast }=(L^{\ast },A^{\ast },A_p^{\ast })$ can be determined by the distribution of the roots of Equation (10(10) $f(\lambda )={\lambda }^3+C_1{\lambda }^2+C_2\lambda +C_3=0,$(10) ). That is, if all the roots of Equation (10(10) $f(\lambda )={\lambda }^3+C_1{\lambda }^2+C_2\lambda +C_3=0,$(10) ) have negative real parts, then $E^{\ast }$ is locally asymptotically stable; if at least one root of Equation (10(10) $f(\lambda )={\lambda }^3+C_1{\lambda }^2+C_2\lambda +C_3=0,$(10) ) has positive real parts, then $E^{\ast }$ is unstable; if any root has zero real part and other roots all have negative real parts, then the stability of $E^{\ast }$ cannot be determined by the linearized system directly.

The following theorem provides a global result on dynamics of the proposed model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) regarding when a colony will collapse.

Theorem 3.2

Extinction of species

If $0<a<{\hat{a}}^{\ast }$ and ${\theta }_0<\frac{{\alpha }_1^3\gamma +4d{\alpha }_1{\beta }^2(1-{\alpha }_1)}{4d{\beta }^2(1-{\alpha }_1{)}^2},$ then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has global stability at $E_0=(0,0,0)$.

Biological implications: Theorem 3.2 has stronger result than results stated in Proposition 3.1 and indicates that the portion of the division of labour invested on larvae a and the nutrient ${\theta }_0$ are important factors determining whether larvae and adult worker ants can survive. This theorem provides a sufficient condition leading to the collapse of the colony.

Theorem 3.3

Stability Conditions

For Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ),

1. Assume that $a>\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}$ and ${\theta }_0<\frac{{\alpha }_1}{1-{\alpha }_1}$, then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has a unique interior equilibrium $E_2=(L_2,A_2,A_{p2})=(\frac{d}{\beta }{A}_2^2,A_2,\frac{{\alpha }_1\beta A_2+{\alpha }_1{A}_2^2}{\beta +{\alpha }_1{A}_2})$. If it satisfies $C_1(E_2)C_2(E_2)>C_3(E_2)>0$, then $E_2$ is locally asymptotically stable.

2. Assume that ${\hat{a}}^{\ast }<a<\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}$ and ${\theta }_0<\frac{{\alpha }_1}{1-{\alpha }_1}$, Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two interior equilibria $E_i=(L_i,A_i,A_{pi})=(\frac{d}{\beta }{A}_i^2,A_i,\frac{{\alpha }_1\beta A_i+{\alpha }_1{A}_i^2}{\beta +{\alpha }_1{A}_i}),i=1,2$, where $E_1<E_2$, if $C_1(E_1)C_2(E_1)-C_3(E_1)<0$ but $C_1(E_2)C_2(E_2)>C_3(E_2)>0$, then the interior equilibrium $E_2$ is locally asymptotically stable while $E_1$ is unstable.

Biological Implications: The results in Lemma 3.1, Theorems 3.1 and 3.3, imply that the division of labour invested on larvae a decreases past the critical point ${\hat{a}}^{\ast }=\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^4{\gamma }^2+4d{\alpha }_1{\beta }^2\gamma (1-{\alpha }_1)({\alpha }_1-(1-{\alpha }_1){\theta }_0)}$ shown in (6(6) $\begin{array}{rl}{\hat{a}}^{\ast }& =\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^4{\gamma }^2+4d{\beta }^2{\alpha }_1\gamma (1-{\alpha }_1)[{\alpha }_1-{\theta }_0(1-{\alpha }_1)]}\\ & =\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}\frac{4d{\beta }^2(1-{\alpha }_1)}{\frac{{\alpha }_1^3\gamma }{{\alpha }_1-(1-{\alpha }_1){\theta }_0}+4d{\beta }^2(1-{\alpha }_1)}<\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}\end{array}$(6) ) and the first dotted line in Figure 2. Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) exhibits a backward bifurcation shown in Figure 2 where $b=0.1,d=0.1,{\alpha }_1=0.3,\beta =0.7,\gamma =0.9$. In Figure 2(a), we set ${\theta }_0=0.1$ and in Figure 2(b), we set ${\theta }_0=0.2$. Based on the expression of the critical value ${\hat{a}}^{\ast }$ shown in (6(6) $\begin{array}{rl}{\hat{a}}^{\ast }& =\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^4{\gamma }^2+4d{\beta }^2{\alpha }_1\gamma (1-{\alpha }_1)[{\alpha }_1-{\theta }_0(1-{\alpha }_1)]}\\ & =\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}\frac{4d{\beta }^2(1-{\alpha }_1)}{\frac{{\alpha }_1^3\gamma }{{\alpha }_1-(1-{\alpha }_1){\theta }_0}+4d{\beta }^2(1-{\alpha }_1)}<\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}\end{array}$(6) ), ${\hat{a}}^{\ast }$ is an increasing function of ${\theta }_0$, which is reflected in the difference between Figure 2(a) and Figure 2(b). The value of ${\theta }_0$ measures the minimum ratio of protein to carbohydrates that can allow the survival of larvae. Simulations shown in Figure 2(a, b) and 3 suggest that the larger value of ${\theta }_0$, the more likely the colony can survive with a larger population of workers A and thus the higher nutrient ratio $\frac{A_p}{A_c}$. In summary, our theoretical work combined with the related simulations suggest that the division of labour invested on larvae a and the minimal nutrient ratio ${\theta }_0$ can affect colony survival, the distribution of the brood and the total forager population affect the protein forager population. For instance, the larger ${\theta }_0$, the more division of labour invested on larvae is required to ensure survival of the colony. Also, under this scenario, the population distribution of brood and workers is smaller. Moreover, Figure 3 shows the bifurcation diagrams of the ratio of $\frac{A_p}{A_c}$ versus the division of labour invested on larvae a with different values of ${\theta }_0$, other parameters values are taken as those in Figure 2.

Figure 2. Backward bifurcation diagram of the division of labour invested on brood a v.s. the total forager population A. Other parameters values are $b=0.1,d=0.1,{\alpha }_1=0.3,\beta =0.7,\gamma =0.9$. The solid line indicates that the equilibrium is locally asymptotically stable while the dashed line indicates that the equilibrium is unstable. The first vertical dotted line is the critical point ${\hat{a}}^{\ast }$ for saddle node bifurcation and the second dotted line is the transition point when the system has two interior equilibrium to one interior equilibrium. The blue colour indicates $E_2$ which is always stable; the green colour indicates $E_1$ which is always unstable; and the red colour is the extinction equilibrium $E_0$. (a) ${\theta }_0=0.1$. (b) ${\theta }_0=0.2$

Figure 3. Bifurcation diagrams of the ratio of $\frac{A_p}{A_c}$ v.s. the division of labour invested on larvae a with: (a) different values of nutrient threshold ${\theta }_0$ and (b) different values of optimal nutrient ratio ${\theta }_m$ when ${\theta }_c=7.8$. Other parameters values are $b=0.1,d=0.1,{\alpha }_1=0.3,\beta =0.7,\gamma =0.9$. (a) Effects on ${\theta }_0$. (b) Effects on ${\theta }_m$.

All our theoretical results can apply to the symmetric case when ${\alpha }_1={\alpha }_2=\alpha$ and ${\theta }_0=2{\theta }_m-{\theta }_c$. Now we focus on the special case of the symmetric case ${\theta }_0=0$. For convenience, let (12) $\begin{array}{r}{a}^{\ast }\triangleq \frac{4b{d}^2{\beta }^2(1-\alpha {)}^2}{{\alpha }^2\gamma [{\alpha }^2\gamma +4d{\beta }^2(1-\alpha )]}=\frac{bd(1-\alpha )}{{\alpha }^2\gamma }\frac{4d{\beta }^2(1-\alpha )}{[{\alpha }^2\gamma +4d{\beta }^2(1-\alpha )]}<\frac{bd(1-\alpha )}{{\alpha }^2\gamma }.\end{array}$(12) Let $\mathrm{\Delta }=a^2{\alpha }^2\gamma [{\alpha }^2\gamma +4d{\beta }^2(1-\alpha )]-4ab{d}^2{\beta }^2(1-\alpha {)}^2$ and $\begin{array}{rl}{A}_1& =\frac{a{\alpha }^2\gamma -\sqrt{\mathrm{\Delta }}}{2a\beta d(1-\alpha )},A_2=\frac{a{\alpha }^2\gamma +\sqrt{\mathrm{\Delta }}}{2a\beta d(1-\alpha )}.\\ (L_{\ast },A_{\ast },{A_p}_{\ast })& =(\frac{d}{\beta }{A}_{\ast }^2,\frac{{\alpha }^2\gamma }{2\beta d(1-\alpha )},\frac{\alpha \beta A_{\ast }+\alpha A_{\ast }^2}{\beta +\alpha A_{\ast }}).\end{array}$ Define $a_1,a_2$ and M as follows: (13) $\begin{array}{rl}{a}_1& =\frac{b({\alpha }^2\gamma +4d{\beta }^2(1-\alpha )-\alpha \sqrt{\gamma ({\alpha }^2\gamma +8d{\beta }^2(1-\alpha ))})}{2{\alpha }^2{\beta }^2\gamma },\\ a_2& =\frac{b({\alpha }^2\gamma +4d{\beta }^2(1-\alpha )+\alpha \sqrt{\gamma ({\alpha }^2\gamma +8d{\beta }^2(1-\alpha ))})}{2{\alpha }^2{\beta }^2\gamma },\\ M& =\frac{1}{{\beta }^2}(2\alpha d{A_2}^2+(\alpha +2d)\beta A_2+3{\beta }^2)+\frac{d}{a\alpha \gamma }(\frac{A_2}{{A_p}_2}(a{A_2}^2-b)+2b).\end{array}$(13)

Theorem 3.4

Dynamics of the Special Symmetric Case

Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is positive invariant in ${\mathbb{R}}_{+}^3$ and every trajectory attracts to a compact set $\mathbb{C}=[0,\frac{\alpha \gamma {\theta }_m}{\beta }]\times [0,\frac{\sqrt{d\alpha \gamma {\theta }_m}}{d}]\times [0,\frac{\alpha \beta \sqrt{d\alpha \gamma {\theta }_m}+{\alpha }^2\gamma {\theta }_m}{d\beta }]$. In addition,

1.	If $a<min\{a^{\ast },\frac{bd}{\alpha \gamma {\theta }_m}\}$, then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has global stability at $E_0=(0,0,0)$.
2.	If $a>\frac{bd(1-\alpha )}{{\alpha }^2\gamma }$, $a_1<a<a_2$ and $\begin{array}{rl}& max\{0,\frac{\beta }{\beta -\alpha A_2}\}\\ & <\frac{A_2}{{A_c}_2}<min\{\frac{\alpha L_2+d{A}_2}{\alpha L_2+d{A_p}_2}[\frac{d(a{A_2}^2-b)}{a\alpha \gamma }+\frac{2{A_p}_2}{A_2}(\frac{bd}{a\alpha \gamma }+1)],M\},\end{array}$ then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has a unique interior equilibrium $E_2=(L_2,A_2,{A_p}_2)$ that is locally asymptotically stable.
3.	If $a^{\ast }<a<\frac{bd(1-\alpha )}{{\alpha }^2\gamma }$, $a_1<a<a_2$, and $\begin{array}{rl}& \frac{A_1}{{A_c}_1}>\frac{1}{{\beta }^2}(2\alpha d{A_1}^2+(\alpha +2d)\beta A_1+3{\beta }^2)-\frac{d}{a\alpha \gamma }(\frac{A_1}{{A_p}_1}(b-a{A_1}^2)-2b),\\ & \frac{A_2}{{A_c}_2}<\frac{1}{{\beta }^2}(2\alpha d{A_2}^2+(\alpha +2d)\beta A_2+3{\beta }^2)+\frac{d}{a\alpha \gamma }(\frac{A_2}{{A_p}_2}(a{A_2}^2-b)+2b),\end{array}$ then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two interior equilibria $E_i=(L_i,A_i,{A_p}_i)=(\frac{d}{\beta }{A}_i^2,A_i,\frac{\alpha \beta A_i+\alpha A_i^2}{\beta +\alpha A_i}),i=1,2$, where the interior equilibrium $E_2$ is locally asymptotically stable while $E_1$ is unstable.
4.	If $a>\frac{bd(1-\alpha )}{{\alpha }^2\gamma }$, $a_1<a<a_2$ and $\begin{array}{rl}& max\{0,\frac{\beta }{\beta -\alpha A},\frac{\beta +\alpha A}{\beta (1-\alpha )}\}\\ & <\frac{A}{A_c}<min\{\frac{\alpha L+dA}{\alpha L+d{A}_p}[\frac{d(a{A}^2-b)}{a\alpha \gamma }+\frac{2A_p}{A}(\frac{bd}{a\alpha \gamma }+1)],\\ & M,\frac{ad{A}^2+bd+a\alpha \gamma }{a\alpha \gamma }\},\end{array}$ then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has a unique interior equilibrium $E_2$ which is globally stable.

Remark 3.3

Theoretical results and numerical simulations (see Section 4) confirm that the special symmetric case of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ), i.e. ${\theta }_0=0$ and ${\theta }_c=2{\theta }_m$, undergoes a backward bifurcation as a decreases past the critical value $a^{\ast }$ defined in (12(12) $\begin{array}{r}{a}^{\ast }\triangleq \frac{4b{d}^2{\beta }^2(1-\alpha {)}^2}{{\alpha }^2\gamma [{\alpha }^2\gamma +4d{\beta }^2(1-\alpha )]}=\frac{bd(1-\alpha )}{{\alpha }^2\gamma }\frac{4d{\beta }^2(1-\alpha )}{[{\alpha }^2\gamma +4d{\beta }^2(1-\alpha )]}<\frac{bd(1-\alpha )}{{\alpha }^2\gamma }.\end{array}$(12) ). Conditions shown in Theorem 3.4 suggest the importance of the protein to carbohydrate ratio, i.e. $\frac{A_2}{{A_c}_2}=1+\frac{{A_p}_2}{{A_c}_2}$, in determining the colony population dynamics.

Summary of Dynamics: According to our analytical results shown in this section, we can conclude that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) undergoes a backward bifurcation as a decreases past ${\hat{a}}^{\ast }$ (or $a^{\ast }$ in the case of ${\theta }_0=0$ and ${\theta }_c=2{\theta }_m$). More specifically, it exhibits the following global dynamics:

1.	If $0<a<{\hat{a}}^{\ast }$ and ${\theta }_0<\frac{{\alpha }_1^3\gamma +4d{\alpha }_1{\beta }^2(1-{\alpha }_1)}{4d{\beta }^2(1-{\alpha }_1{)}^2}$, then the colony collapses due to the lack of efforts of division of labour invested on larvae and the minimum nutrient requirement ${\theta }_0$ being too low.
2.	If $0<{\hat{a}}^{\ast }<a<\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}$, the survival of the colony depends on its initial population size.
3.	If $a>\frac{bd(1-{\alpha }_1)}{{\alpha }_1\gamma ({\alpha }_1-(1-{\alpha }_1){\theta }_0)}>0$, then the colony persists.

Our theoretical results suggest that the survival rate of larva to worker $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$ plays critical roles in determining colony population dynamics. We assume that $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$ takes the form of (1(1) $S_{L_{max}}\left(\frac{A_p}{A_c}\right)=\{\begin{array}{ll}0& \text{ when }0\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_0,\\ {\alpha }_1({\displaystyle \frac{A_p}{A_c}}-{\theta }_0)& \text{when }{\theta }_0\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_m,\\ -{\alpha }_2({\displaystyle \frac{A_p}{A_c}}-{\theta }_c)& \text{when }{\theta }_m\le {\displaystyle \frac{A_p}{A_c}}<{\theta }_c,\\ 0& \text{when }{\theta }_c\le {\displaystyle \frac{A_p}{A_c}},\end{array}$(1) ) based on relevant biological studies. Our analysis implies that the values of parameters ${\alpha }_1$ and ${\theta }_0$ in $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$ have pronounced impacts on dynamical outcomes. In the next section, we use bifurcation diagrams to explore detailed impacts.

4. Numerical simulations

In this section, we use numerical simulations to illustrate equilibrium dynamics of the proposed model and obtain further biological insights on the dynamical outcomes of certain life history parameters of the colony.

For the general case of $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$, the dynamics of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) depends on the division of labour a, egg laying rate γ, the scaling factor on the brood survival rate due to the nutritional status ${\alpha }_1$, the minimal nutrition ratio ${\theta }_0$, the maturation rate β, and the natural mortality d. To explore the effects of a and ${\theta }_0$, we perform bifurcation diagrams in Figures 2 and  3 by setting $b=0.1,d=0.1,{\alpha }_1=0.3,\beta =0.7,\gamma =\mathrm{0.9.}$Figures 2 and  3 suggest that (1) small values of division of labour a can lead to colony collapse; (2) intermediate values of a can make the system go through saddle node bifurcation; and (3) large values of a can insure colony survival. This implies that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) goes through backward bifurcation on a. We can see that the larger value of a can lead to the larger population A (see Figure 2) and the larger nutrient ratio $\frac{A_p}{A_c}$ (see Figure 3). Figures 2 and 3 also show the effects of the minimal nutrient requirement for brood survival ${\theta }_0$: The larger value of ${\theta }_0$, (1) the larger critical threshold ${\hat{a}}^{\ast }$; (2) the smaller population A; and (3) the smaller nutrient status, i.e. the smaller value of $\frac{A_p}{A_c}$.

Next, we perform bifurcation diagrams of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) regarding how the minimum nutritional requirement ${\theta }_0$ and the scaling factor of survival probability of brood ${\alpha }_1$ affect population dynamics of the colony in Figure 4. Figure 4(a) shows that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) exhibits reversed backward bifurcation on ${\theta }_0$. Figure 4(b) suggests that the larger value of ${\alpha }_1$, the larger population of worker A and the better probability of colony survival.

Figure 4. (a) The bifurcation diagrams of $L,A$ and $A_p$ v.s. the nutrient threshold ${\theta }_0$ with ${\alpha }_1=0.3$; and (b) the bifurcation diagram of A v.s. the nutrient threshold ${\theta }_0$ with different values of ${\alpha }_1$. Other parameters values are taken as $a=0.15,b=0.1,d=0.1,\beta =0.7,\gamma =0.9$. (a) Effects of ${\theta }_0$ on population dynamics. (b) Effects of ${\alpha }_1$ on A

In the remaining of this section, we focus on the symmetric case of $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$ shown in (3(3) $\begin{array}{rl}& \frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}=\alpha ;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}=-\alpha {\left(\frac{A_p}{A_c}\right)}^2\text{ ,when, }2{\theta }_m-{\theta }_c\le \frac{A_p}{A_c}<{\theta }_m,\\ & \frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}=-\alpha ;\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}=\alpha {\left(\frac{A_p}{A_c}\right)}^2\text{ ,when, }{\theta }_m\le \frac{A_p}{A_c}<{\theta }_c.\end{array}$(3) ) where ${\alpha }_1={\alpha }_2=\alpha$ and ${\theta }_0=2{\theta }_m-{\theta }_c$.

Special symmetric case ${\theta }_c=2{\theta }_m$ (i.e. ${\theta }_0=0$): Figure 5(a) provides an example of bifurcation diagram on division of labour invested on larvae a of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) by choosing the following parameters values: $b=0.1,d=0.1,\alpha =0.3,\beta =0.7,\gamma =0.9,{\theta }_c=8,{\theta }_m=4.$ Figure 5(a) shows that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) goes through backward bifurcation at $a^{\ast }=\frac{4b{d}^2{\beta }^2(1-\alpha {)}^2}{{\alpha }^2\gamma [{\alpha }^2\gamma +4d{\beta }^2(1-\alpha )]}=0.054$. If $a<a^{\ast }$, then colony collapses; if $0.054=a^{\ast }<a<\frac{bd(1-\alpha )}{{\alpha }^2\gamma }$, Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two positive interior equilibria $E_1=(L_1,A_1,{A_p}_1)$ and $E_2=(L_2,A_2,{A_p}_2)$ with $E_2$ being locally stable; and if $a>\frac{bd(1-\alpha )}{{\alpha }^2\gamma }$, then the colony survives.

Figure 5. Bifurcation diagrams of the division of labour invested on larvae a for Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ). Backward bifurcation occurs at $a^{\ast }\triangleq \frac{4b{d}^2{\beta }^2(1-\alpha {)}^2}{{\alpha }^2\gamma [{\alpha }^2\gamma +4d{\beta }^2(1-\alpha )]}$ for the case of $2{\theta }_m={\theta }_c=8$ (a); and at ${\tilde{a}}^{\ast }=\frac{4b{d}^2{\beta }^2(1-\alpha {)}^2}{{\alpha }^4{\gamma }^2+4d\alpha {\beta }^2\gamma (1-\alpha )(\alpha -(1-\alpha )(2{\theta }_m-{\theta }_c))}$ for the case of $2{\theta }_m>{\theta }_c=7.8$ (b). Other parameters values are taken as $b=0.1,d=0.1,\alpha =0.3,\beta =0.7,\gamma =0.9,{\theta }_m=4$. The solid line indicates that the equilibrium is locally asymptotically stable while the dotted line indicates that the equilibrium is unstable. The blue colour indicates $E_2$ which is always stable; the green colour indicates $E_1$ which is always unstable; and the red colour is the extinction equilibrium $E_0$. (a) $2{\theta }_m={\theta }_c=8$. (b) $2{\theta }_m>{\theta }_c=7.8$

Figure 6 provides two examples of population dynamics of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) to show the effects of a. The initial condition is $(L(0),A(0),A_p(0))=(0.09,0.1,0.035)$ and other parameters values are the same as in Figure 5(a). According to Theorem 3.4, Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has global stability at $E_0=(0,0,0)$ (i.e. colony collapses) if $a=0.07<min\{a^{\ast },\frac{bd}{\alpha \gamma {\theta }_m}\}$ (see Figure 6(a)) and Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two positive interior equilibria with $E_2=(0.434,1.743,1.045)$ being a locally asymptotically stable interior equilibrium if $a=0.1>a^{\ast }$ (see Figure 6(b)). It indicates that larvae L, worker ants A, and the ants collecting proteinaceous material $A_p$ can coexist with proper initial conditions if a is in the intermediate range.

Figure 6. The time series of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) when ${\theta }_c=2{\theta }_m=8$ with an initial value $(L(0),A(0),A_p(0))=(0.09,0.1,0.035)$ and other parametric values being $b=0.1,d=0.1,\alpha =0.3,\beta =0.7,\gamma =0.9,$ which is the same set of parameters values in Figure 5(a). Figure (a) is the case when a = 0.07 where population goes extinct over time. Figure (b) is the case when a = 0.1 where the colony has a locally asymptotically stable interior equilibrium $E_2=(0.434,1.743,1.045)$. (a) Time series of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) when ${\theta }_c=2{\theta }_m$ when the portion of the division of labour invested on larvae is a = 0.07. (b) The time series of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) when ${\theta }_c=2{\theta }_m$ when the portion of the division of labour invested on larvae is a = 0.1.

Symmetrical case ${\theta }_c\le 2{\theta }_m$ (i.e. ${\theta }_0=2{\theta }_m-{\theta }_c>0$): Figure 5(b) provides an example of bifurcation diagram on division of labour invested on larvae (a) of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) by taking same parameter values in Figure 5(a) except that ${\theta }_c=7.8$. Figure 5(b) shows that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) undergoes a bifurcation as the portion of the division of labour invested on larvae a decreasing past ${\tilde{a}}^{\ast }=\frac{4b{d}^2{\beta }^2(1-\alpha {)}^2}{{\alpha }^4{\gamma }^2+4d\alpha {\beta }^2\gamma (1-\alpha )(\alpha -(1-\alpha )(2{\theta }_m-{\theta }_c))}$. Notice that the symmetrical case has ${\theta }_0=2{\theta }_m-{\theta }_c>0$. Thus, the comparisons between Figure 5(a) and Figure 5(b) can provide insights on the effect of a and ${\theta }_c$ (or ${\theta }_0$ because we set ${\theta }_m=4$ and ${\theta }_0=2{\theta }_m-{\theta }_c$): (1) The smaller value of ${\theta }_c$, the larger critical threshold $a^{\ast }$; (2) The smaller value of ${\theta }_c$, the smaller population A. The dynamical outcomes of symmetrical cases are similar to the general case shown in Figures 2 and  3.

To understand the effects of the optimal nutrient ratio ${\theta }_m$ (or ${\theta }_0=2{\theta }_m-{\theta }_c=2{\theta }_m-7.8$), we perform a bifurcation diagram on ${\theta }_m$ shown in Figure 7 by setting $a=0.15,b=0.1,d=0.1,\alpha =0.3,\beta =0.7,\gamma =0.9,{\theta }_c=\mathrm{7.8.}$ Notice that ${\theta }_0=2{\theta }_m-{\theta }_c\ge 0$, thus the value of ${\theta }_m$ in Figure 7 starts with ${\theta }_m=1/2{\theta }_c=3.9$. Figure 7 shows that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) exhibits reversed backward bifurcation on ${\theta }_m$ (or ${\theta }_0$): (1) small values of the optimal nutrient ratio ${\theta }_m$ can insecure the persistence of the colony; (2) intermediate values of ${\theta }_m$ can go through saddle node bifurcation; and (3) large values of ${\theta }_m$ can lead to colony collapse.

Figure 7. Bifurcation diagram for Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) on the optimal nutrient ratio ${\theta }_m$. An unstable interior equilibrium (green dotted) bifurcates ${\theta }_m^{\ast }=\frac{1}{2}({\theta }_c+\frac{\alpha }{1-\alpha }+\frac{{\alpha }^3\gamma }{4d{\beta }^2(1-\alpha {)}^2}-\frac{bd}{a\alpha \gamma })$. The solid line indicates that the equilibrium is stable, while the dotted line indicates that the equilibrium is unstable. Parameters values: $a=0.15,b=0.1,d=0.1,\alpha =0.3,\beta =0.7,\gamma =0.9,{\theta }_c=7.8$.

5. Conclusion

Variation in nutrient consumption among individuals is considered a conserved mechanism regulating castes and division of labour in social insects colonies. In eusocial insects, foragers, who perform food collection tasks, need to satisfy their own nutrient requirements in addition to those of the non-foraging workers, as well as the larvae and queen(s), which have significantly higher protein needs [14]. In this paper, we propose and study a nonlinear differential equations system to explore how nutritional status may regulate population dynamics and foraging task allocation of social insect colonies by applying adaptive modeling framework. Our model assumes that foragers adjust their preferences in favour of food sources containing limiting nutrients to maintain colony growth and reproduction [7,11,26].

Our proposed model consists of a population of larvae L, foragers collecting carbohydrate $A_c$, and foragers collecting protein $A_p$. We assume that the survival rate of larvae is determined by the available nutrition in the colony, which is reflected through the ratio of workers collecting protein to those collecting carbohydrates $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$. Our formulation of $S_{L_{max}}\left(\frac{A_p}{A_c}\right)$ is based on biological studies (see  [11,26]) and embeds with an adaptive modeling approach adopted from [15,16,19]. Our theoretical results and bifurcation analysis conclude that our proposed model exhibits backward bifurcations that generate bistability (see examples in Figure 6). The bistability of the colony implies that initial conditions are important for colony survival under certain ranges of life history parameters. More specifically, the dynamical features and the related biological implications of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) can be summarized as follows:

1. The nutrition status measured by $\frac{A_p}{A_c}$ is an increasing function of the total population of workers A, or vice versa. This result may stem from the assumption that larvae (or brood in general) have higher nutritional needs to ensure survival. The biological implications of this are that higher nutritional status $\frac{A_p}{A_c}$ can lead to a better survival rate of larvae, thus the colony can grow with larger worker population A.

2. The survival probability of brood is an increasing function of the following important life history parameters of colony:

   1. The division of labour invested on brood measured by a can have huge impacts on dynamical outcomes of the colony. From Lemma 3.1, Theorems 3.1 and 3.3, Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) exhibits a backward bifurcation (shown in Figure 2) as a decreasing past its critical point ${\hat{a}}^{\ast }=\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^4{\gamma }^2+4d{\alpha }_1{\beta }^2\gamma (1-{\alpha }_1)({\alpha }_1-(1-{\alpha }_1){\theta }_0)}$ which is an increasing function of the maturation rate β, the minimal nutrient ratio ${\theta }_0$; and a decreasing function of queen(s) laying egg rate γ. The larger the value of a, the more likely it is that the colony survives and grows.

   2. Effects of nutrient thresholds ${\theta }_0,{\theta }_m,{\theta }_c$: Figures 2, 3 and 7 suggest that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) exhibits reversed backward bifurcation on ${\theta }_m$ (or ${\theta }_0$): (1) small values of the optimal nutrient ratio ${\theta }_m$ can insecure the persistence of the colony; (2) intermediate values of ${\theta }_m$ can go through saddle node bifurcation that leads to bistability; and (3) large values of ${\theta }_m$ can lead to colony collapsing.

      In addition, the larger value of the optimal nutrient ratio ${\theta }_m$ (or ${\theta }_0$) leads to (1) larger critical threshold ${\hat{a}}^{\ast }$; (2) smaller population A; and (3) smaller nutrient status, i.e. smaller value of $\frac{A_p}{A_c}$.

   3. Effects of the brood survival rate ${\alpha }_1$: Figure 4(a) shows that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) exhibits reversed backward bifurcation on ${\theta }_0$. Figure 4(b) suggests that the larger value of ${\alpha }_1$, the larger the population of workers A will be, which increases the probability of colony survival.

Our proposed model and study provide new insights into the strategies used by social insects (such as harvesting ants) facing nutritional challenges, and our results deepen our understanding of their nutritional ecology. Task allocation has been studied in social insects, that is, how colonies change the allocation of tasks in response to changing colony needs. One of future directions would be extending our current model to include more tasks such as brood care, foraging, and study how different tasks are related to colony needs including nutritional requirement. An another future direction is to extend our current model to include an additional level such as food resource of the colony. For example, leaf-cutter ants collect leaves as food resource to cultivate fungi and harvest the fruits of fungi as their food. The nutrient requirement of the leaf-cutter ants colony has two levels: one is the needs of the colony itself such as brood and the other one is the needs of the fungi. It would be interesting to explore how nutrient needs of the colony and fungus garden affect the foraging behaviour of leaf-cutter ants during its ontology.

6. Proofs

Proof of Lemma 3.1

Proof.

For any $L,A\in {\mathbb{R}}_{+}^2$, from Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ), we obtain $L^{\prime }{|}_{L(0)=0}\ge 0$ and $A^{\prime }{|}_{A(0)=0}\ge 0$ for all $t\ge 0$. Since $A=A_p+A_c$, then $A_p^{\prime }{|}_{A_p(0)=0}\ge 0$ for all $t\ge 0$. Moreover, if $L(0)=0,A(0)=0$ and $A_p(0)=0$, then $(L(t),A(t),A_p(t))=(0,0,0)$ for all $t\ge 0$. If $L(0)>0,A(0)>0$ and $A_p(0)>0$, then by continuity arguments, it is impossible for either $L(t)$ or $A(t)$ or $A_p(t)$ to drop below 0. Hence, for any $L(0)\ge 0,A(0)\ge 0$ and $A_p(0)\ge 0$, we obtain $L(t)\ge 0,A(t)\ge 0$ and $A_p(t)\ge 0$ for all $t\ge 0$.

Now assume $L(0)\ge 0,A(0)\ge 0$ and $A_p(0)\ge 0$, then since the function of $S_L\left(\frac{A_p}{A_c}\right)$ exists maximum when $\frac{A_p}{A_c}\in (max\{0,{\theta }_0\},{\theta }_c)$ and according to the expression of $L^{\prime }$, we have $L^{\prime }={\alpha }_1\gamma (\frac{A_p}{A_c}-{\theta }_0)\frac{a{A}^2}{b+a{A}^2}-\beta L\le {\alpha }_1\gamma ({\theta }_m-{\theta }_0)-\beta L$ for all $t\ge 0$ when ${\theta }_m>{\theta }_0$. Thus, a standard comparison theorem shows that $\underset{t\to \mathrm{\infty }}{lim sup}L(t)\le \frac{{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}{\beta }$. This indicates that for any $ϵ>0$, there exists T large enough, such that $L(t)\le \frac{{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}{\beta }+ϵ\text{for all}t>T.$ Therefore, from the expression of $A^{\prime }$, we have $A^{\prime }=\beta L-d{A}^2\le \beta (\frac{{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}{\beta }+ϵ)-d{A}^2\text{for all}t>T.$ Since ε can be arbitrarily small, thus $\underset{t\to \mathrm{\infty }}{lim sup}A(t)\le \frac{\sqrt{d{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}}{d}$. Thus, we have shown that the Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is positively invariant and bounded in ${\mathbb{R}}_{+}^2$. More specifically, the compact set $[0,\frac{{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}{\beta }]\times [0,\frac{\sqrt{d{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}}{d}]$ attracts all points in ${\mathbb{R}}_{+}^2$. Due to $A=A_p+A_c$ and the boundedness of $A(t)$, hence we can obtain $A_p$ is bounded for all $t\ge 0$.

Moreover, if $L(0)>0,A(0)>0$ and $A_p(0)>0$, then we have follows: $\begin{array}{rl}{L}^{\prime }& ={\alpha }_1\gamma (\frac{A_p}{A_c}-{\theta }_0)\frac{a{A}^2}{b+a{A}^2}-\beta L\ge -\beta L\Rightarrow L(t)\ge L(0){\mathrm{e}}^{-\beta t}>0,\\ A^{\prime }& =\beta L-d{A}^2\ge -d{A}^2\Rightarrow A(t)\ge \frac{A(0)}{1+\mathrm{d}t}>0,\\ A_p^{\prime }& ={\alpha }_1\beta L+{\alpha }_1{A}_{c}L-d{A}_{p}A\ge -d{A}_{p}A\Rightarrow A_p(t)\ge A_p(0)e^{-d{\int }_0^{t}A(s)\mathrm{d}s}>0.\end{array}$ Therefore, if $L(0)>0,A(0)>0$ and $A_p(0)>0$, then $L(t)>0,A(t)>0$ and $A_p(t)>0$ for all t>0.

Proof of Theorem 3.1

Proof.

It is easy to see that $E_0=(0,0,0)$ is always an equilibrium of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ). The nullclines of (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) can be found as $\begin{array}{rl}{L}^{\prime }& =0⟹{\alpha }_1\gamma (\frac{A_p}{A_c}-{\theta }_0)\frac{a{A}^2}{b+a{A}^2}-\beta L=0,\\ A^{\prime }& =0⟹L=\frac{d}{\beta }{A}^2,\\ A_p^{\prime }& =0⟹{\alpha }_1\beta L+{\alpha }_1{A}_{c}L-d{A}_{p}A=0⟹L=\frac{d{A}_{p}A}{{\alpha }_1\beta +{\alpha }_1{A}_c}.\end{array}$ By solving $\frac{d{A}_{p}A}{{\alpha }_1\beta +{\alpha }_1{A}_c}=\frac{d}{\beta }{A}^2$ for $A_p$, we have $A_p=\frac{{\alpha }_1\beta A+{\alpha }_1{A}^2}{\beta +{\alpha }_{1}A}$ and substitute it to $L^{\prime }=0$, which results in the following equation: (14) $ad\beta (1-{\alpha }_1)A^2-a{\alpha }_1^2\gamma A+\beta [(1-{\alpha }_1)(bd+a{\alpha }_1\gamma {\theta }_0)-a{\alpha }_1^2\gamma ]=0.$(14) The roots of (14(14) $ad\beta (1-{\alpha }_1)A^2-a{\alpha }_1^2\gamma A+\beta [(1-{\alpha }_1)(bd+a{\alpha }_1\gamma {\theta }_0)-a{\alpha }_1^2\gamma ]=0.$(14) ) are given by $A_1=\frac{a{\alpha }_1^2\gamma -\sqrt{\mathrm{\Delta }}}{2a\beta d(1-{\alpha }_1)},A_2=\frac{a{\alpha }_1^2\gamma +\sqrt{\mathrm{\Delta }}}{2a\beta d(1-{\alpha }_1)},$ where $\mathrm{\Delta }=a(a{\alpha }_1^4{\gamma }^2-4d{\beta }^2[(1-{\alpha }_1{)}^2(bd+a{\alpha }_1\gamma {\theta }_0)-a{\alpha }_1^2\gamma (1-{\alpha }_1)])$.

Thus, we have the following three cases:

Let ${\theta }_0^{\ast }=\frac{{\alpha }_1^3\gamma }{4d{\beta }^2(1-{\alpha }_1{)}^2}+\frac{{\alpha }_1}{1-{\alpha }_1}-\frac{bd}{a{\alpha }_1\gamma }$.

1. If ${\theta }_0>{\theta }_0^{\ast }$ and $a>\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^2\gamma ({\alpha }_1^2\gamma +4d{\beta }^2(1-{\alpha }_1))}$, then there is only one trivial equilibrium: $E_0=(0,0,0)$ and no other positive interior equilibrium.

2. If ${\theta }_0={\theta }_0^{\ast }$ and $a>\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^2\gamma ({\alpha }_1^2\gamma +4d{\beta }^2(1-{\alpha }_1))}$, then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two positive equilibria which collapse into one equilibrium $E_{\ast }$ as $(L_{\ast },A_{\ast },{A_p}_{\ast })=(\frac{d}{\beta }{A}_{\ast }^2,\frac{{\alpha }_1^2\gamma }{2\beta d(1-{\alpha }_1)},\frac{{\alpha }_1\beta A_{\ast }+{\alpha }_1{A}_{\ast }^2}{\beta +{\alpha }_1{A}_{\ast }}).$ Or if $0<{\theta }_0<\frac{{\alpha }_1}{1-{\alpha }_1}-\frac{bd}{a{\alpha }_1\gamma }$ and $a>\frac{bd(1-{\alpha }_1)}{{\alpha }_1^2\gamma }$, then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has only one positive equilibrium $(L_2,A_2,A_{p2})=(\frac{d}{\beta }{A}_2^2,A_2,\frac{{\alpha }_1\beta A_2+{\alpha }_1{A}_2^2}{\beta +{\alpha }_1{A}_2}).$

3. If $max\{0,\frac{{\alpha }_1}{1-{\alpha }_1}-\frac{bd}{a{\alpha }_1\gamma }\}<{\theta }_0<{\theta }_0^{\ast }$ and $a>\frac{bd(1-{\alpha }_1)}{{\alpha }_1^2\gamma }$, then Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two positive equilibria in the following form: $\begin{array}{rl}(L_1,A_1,A_{p1})& =(\frac{d}{\beta }{A}_1^2,A_1,\frac{{\alpha }_1\beta A_1+{\alpha }_1{A}_1^2}{\beta +{\alpha }_1{A}_1})\mathrm{and}\\ (L_2,A_2,A_{p2})& =(\frac{d}{\beta }{A}_2^2,A_2,\frac{{\alpha }_1\beta A_2+{\alpha }_1{A}_2^2}{\beta +{\alpha }_1{A}_2}).\end{array}$

Proof of Theorem 3.2

Proof.

From Proposition 3.1, we know that for some initial condition taken in ${\mathbb{R}}_{+}^3$ and if $a<\frac{bd}{{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}$, the trajectory of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is converging to the origin $E_0=(0,0,0)$. And according to Theorem 3.1, if $0<a<{\hat{a}}^{\ast }$ and ${\theta }_0<\frac{{\alpha }_1^3\gamma +4d{\alpha }_1{\beta }^2(1-{\alpha }_1)}{4d{\beta }^2(1-{\alpha }_1{)}^2}$, then there is only one trivial equilibrium $E_0=(0,0,0)$ and no other positive equilibrium. Therefore, we can conclude that Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has global stability at $(0,0,0)$ when $a<min\{{\hat{a}}^{\ast },\frac{bd}{{\alpha }_1\gamma ({\theta }_m-{\theta }_0)}\}$ and ${\theta }_0<\frac{{\alpha }_1^3\gamma +4d{\alpha }_1{\beta }^2(1-{\alpha }_1)}{4d{\beta }^2(1-{\alpha }_1{)}^2}$.

Proof of Theorem 3.3

Proof.

The local stability of equilibria is determined by computing the eigenvalues of the Jacobian matrix about each equilibrium.

Let $E^{\ast }=(L^{\ast },A^{\ast },A_p^{\ast })$ be an arbitrary positive equilibrium of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ). The Jacobian matrix at this equilibrium is (15) $J{|}_{E^{\ast }}=\left(\begin{array}{ccc}-\beta & J_{12}& J_{13}\\ \beta & -2d{A}^{\ast }& 0\\ J_{31}& J_{32}& J_{33}\end{array}\right),$(15) where $\begin{array}{rl}{J}_{12}& =\frac{a{\alpha }_1\gamma A^{\ast }[A_p^{\ast }(b{A}^{\ast }-a{A}^{\ast 3}-2b{A}_p^{\ast })-2b(A^{\ast }-A_p^{\ast }{)}^2{\theta }_0]}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A}^{\ast 2}{)}^2},\\ J_{13}& =\frac{a{\alpha }_1\gamma A^{\ast 3}}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A}^{\ast 2})}>0,\\ J_{31}& ={\alpha }_1(\beta +A^{\ast }-A_p^{\ast })>0,J_{32}={\alpha }_1{L}^{\ast }-d{A}_p^{\ast },J_{33}=-({\alpha }_1{L}^{\ast }+d{A}^{\ast })<0.\end{array}$ Then we have the characteristic equation of $J{|}_{E^{\ast }}$ is (16) $\hat{f}(\lambda )={\lambda }^3+C_1{\lambda }^2+C_2\lambda +C_3=0,$(16) where $\begin{array}{rl}{C}_1& =\beta +\alpha L^{\ast }+3d{A}^{\ast }>0,\\ C_2& =J_{11}{J}_{33}+J_{11}{J}_{22}+J_{22}{J}_{33}-J_{21}{J}_{12}-J_{31}{J}_{13}\\ & =2d\beta A^{\ast }+(\beta +2d{A}^{\ast })({\alpha }_1{L}^{\ast }+d{A}^{\ast })-\frac{a{\alpha }_1^2\gamma A^{\ast 3}(\beta +A^{\ast }-A_p^{\ast })}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A}^{\ast 2})}\\ & -\frac{a{\alpha }_1\beta \gamma A^{\ast }[A_p^{\ast }(b{A}^{\ast }-a{A}^{\ast 3}-2b{A}_p^{\ast })-2b(A^{\ast }-A_P^{\ast }{)}^2{\theta }_0]}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A}^{\ast 2}{)}^2},\\ C_3& =-\det (J{|}_{E_i^{\ast }})=J_{11}{J}_{22}{J}_{33}+J_{21}{J}_{32}{J}_{13}-J_{21}{J}_{12}{J}_{33}-J_{31}{J}_{22}{J}_{13}\\ & =-2d\beta A^{\ast }({\alpha }_1{L}^{\ast }+d{A}^{\ast })\\ & +\frac{a{\alpha }_1\gamma A^{\ast 3}[\beta ({\alpha }_1{L}^{\ast }-d{A}_p^{\ast }+2d{\alpha }_1{A}^{\ast })+2d{\alpha }_1{A}^{\ast }(A^{\ast }-A_p^{\ast })]}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A}^{\ast 2})}\\ & +\frac{a{\alpha }_1\beta \gamma A^{\ast }({\alpha }_1{L}^{\ast }+d{A}^{\ast })[A_p^{\ast }(b{A}^{\ast }-a{A}^{\ast 3}-2b{A}_p^{\ast }-2b(A^{\ast }-A_p^{\ast }{)}^2{\theta }_0)]}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A}^{\ast 2}{)}^2}.\end{array}$ This indicates the following two cases:

1. If Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has a unique interior equilibrium $E_2=(L_2,A_2,A_{p2})=(\frac{d}{\beta }{A}_2^2,A_2,\frac{{\alpha }_1\beta A_2+{\alpha }_1{A}_2^2}{\beta +{\alpha }_1{A}_2})$, then under the conditions $0<{\theta }_0<\frac{{\alpha }_1}{1-{\alpha }_1}-\frac{bd}{a{\alpha }_1\gamma }$ and $a>\frac{bd(1-{\alpha }_1)}{{\alpha }_1^2\gamma }$ and $C_1(E_2)C_2(E_2)>C_3(E_2)>0$, thus, by applying the Routh-Hurwitz criterion, we can obtain that the interior equilibrium $E_2$ of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is locally asymptotically stable.

2. If Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two interior equilibria $E_i=(L_i,A_i,A_{pi})=(\frac{d}{\beta }{A}_i^2,A_i,\frac{{\alpha }_1\beta A_i+{\alpha }_1{A}_i^2}{\beta +{\alpha }_1{A}_i}),i=1,2$ where $E_1<E_2$, then under the conditions $max\{0,\frac{{\alpha }_1}{1-{\alpha }_1}-\frac{bd}{a{\alpha }_1\gamma }\}<{\theta }_0<{\theta }_0^{\ast }$, $a>\frac{4b{d}^2{\beta }^2(1-{\alpha }_1{)}^2}{{\alpha }_1^2\gamma ({\alpha }_1^2\gamma +4d{\beta }^2(1-{\alpha }_1))}$, and $C_1(E_1)C_2(E_1)-C_3(E_1)<0$ but $C_1(E_2)C_2(E_2)>C_3(E_2)>0$, we can obtain that the interior equilibrium $E_2$ is locally asymptotically stable while $E_1$ is unstable.

Proof of Theorem 3.4

Proof.

For any $L,A,A_p\in {\mathbb{R}}_{+}^3$, note that $\begin{array}{rl}{L}^{\prime }{|}_{L=0}& =\alpha \gamma \frac{A_p}{A-A_p}\frac{a{A}^2}{b+a{A}^2}\ge 0,\\ A^{\prime }{|}_{A=0}& =\beta L\ge 0,\\ A_p^{\prime }{|}_{A_p=0}& =\alpha \beta L+\alpha AL\ge 0,\end{array}$ thus according to Theorem A.4 (p. 423) of [31], we can conclude that the model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) is positive invariant in ${\mathbb{R}}_{+}^3$. Now we can proceed to show the boundedness of the system. First, assume $L(0)\ge 0,A(0)\ge 0$ and $A_p(0)\ge 0$, then since the function of $S_L\left(\frac{A_p}{A_c}\right)$ exists maximum when $0<\frac{A_p}{A_c}\le {\theta }_m$ and according to the expression of $L^{\prime }$, we have the following inequalities due to the property of positive invariance: $L^{\prime }=\alpha \gamma \frac{A_p}{A_c}\frac{a{A}^2}{b+a{A}^2}-\beta L\le \alpha \gamma {\theta }_m-\beta L$ which implies that $\underset{t\to \mathrm{\infty }}{lim sup}L(t)\le \frac{\alpha \gamma {\theta }_m}{\beta }.$ This suggests that there exists $ϵ>0$ such that the following inequalities hold as time t is large enough, $A^{\prime }=\beta L-d{A}^2\le \beta (\frac{\alpha \gamma {\theta }_m}{\beta }+ϵ)-d{A}^2$ which indicates that $\underset{t\to \mathrm{\infty }}{lim sup}A(t)\le \sqrt{\frac{\alpha \gamma {\theta }_m}{d}}.$ Then, we also have the following inequalities hold as time t is large enough, $\begin{array}{rl}{A}_p^{\prime }& =\alpha \beta L+\alpha (A-A_p)L-dA{A}_p\le \alpha (\frac{\alpha \gamma {\theta }_m}{\beta }+ϵ)(\beta +\sqrt{\frac{\alpha \gamma {\theta }_m}{d}}+ϵ)\\ & -d(\sqrt{\frac{\alpha \gamma {\theta }_m}{d}}+ϵ)A_p\end{array}$ which shows that $\underset{t\to \mathrm{\infty }}{lim sup}{A}_p(t)\le \alpha (\sqrt{\frac{\alpha \gamma {\theta }_m}{d}}+\frac{\alpha \gamma {\theta }_m}{d\beta }).$ Therefore, every trajectory starting from ${\mathbb{R}}_{+}^3$ converges to the compact set $\mathbb{C}=[0,\frac{\alpha \gamma {\theta }_m}{\beta }]\times [0,\sqrt{\frac{\alpha \gamma {\theta }_m}{d}}]\times [0,\alpha (\sqrt{\frac{\alpha \gamma {\theta }_m}{d}}+\frac{\alpha \gamma {\theta }_m}{d\beta })].$ Let $E^{\ast }=(L^{\ast },A^{\ast },A_p^{\ast })$ be an interior equilibrium of Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ). Then its stability is determined by the eigenvalues ${\lambda }_i(E^{\ast }),i=1,2,3$ of its associated Jacobian matrix as follows: $J{|}_{E^{\ast }}=\left(\begin{array}{ccc}-\beta & -{\displaystyle \frac{a\alpha \gamma A^{\ast }{A}_p^{\ast }(a{A^{\ast }}^3-b{A}^{\ast }+2b{A}_p^{\ast })}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A^{\ast }}^2{)}^2}}& {\displaystyle \frac{a\alpha \gamma {A^{\ast }}^3}{(A^{\ast }-A_p^{\ast }{)}^2(b+a{A^{\ast }}^2)}}\\ \beta & -2d{A}^{\ast }& 0\\ \alpha (\beta +A^{\ast }-A_p^{\ast })& \alpha L^{\ast }-d{A}_p^{\ast }& -\alpha L^{\ast }-d{A}^{\ast }\end{array}\right),$ since $\beta L^{\ast }=d{A^{\ast }}^2,\beta L^{\ast }=\alpha \gamma \frac{A_p^{\ast }}{A_c^{\ast }}\frac{a{A^{\ast }}^2}{b+a{A^{\ast }}^2}$, and $\frac{A^{\ast }}{A_p^{\ast }}=\frac{\beta }{\alpha (\beta +A_c^{\ast })}$. Therefore, we have $J{|}_{E^{\ast }}=\left(\begin{array}{ccc}-\beta & -{\displaystyle \frac{d^2{A^{\ast }}^2(a{A^{\ast }}^2-b)}{a\alpha \gamma A_p^{\ast }}}-\frac{2b{d}^2{A}^{\ast }}{a\alpha \gamma }& {\displaystyle \frac{d{A^{\ast }}^3}{A_p^{\ast }{A}_c^{\ast }}}\\ \beta & -2d{A}^{\ast }& 0\\ \alpha (\beta +A^{\ast }-A_p^{\ast })& \alpha L^{\ast }-d{A}_p^{\ast }& -\alpha L^{\ast }-d{A}^{\ast }\end{array}\right),$ and the characteristic equation of $J{|}_{E^{\ast }}$ is $f(\lambda )={\lambda }^3+c_1{\lambda }^2+c_2\lambda +c_3=0,$ where $\begin{array}{rl}{c}_1& =-\mathrm{tr}(J{|}_{E^{\ast }})=-({\lambda }_1(E^{\ast })+{\lambda }_2(E^{\ast })+{\lambda }_3(E^{\ast }))=\beta +3d{A}^{\ast }+\alpha L^{\ast }>0,\\ c_2& =(\beta +2d{A}^{\ast })(\alpha L^{\ast }+d{A}^{\ast })+\frac{\beta d^2{A}^{\ast }}{a\alpha \gamma }(\frac{A^{\ast }(a{A^{\ast }}^2-b)}{A_p^{\ast }}+2b)+2d\beta A^{\ast }-\frac{\beta d{A^{\ast }}^2}{A_c^{\ast }},\\ c_3& =-\det (J{|}_{E^{\ast }})=-{\lambda }_1(E^{\ast }){\lambda }_2(E^{\ast }){\lambda }_3(E^{\ast })\\ & =[\frac{\beta d^2{A}^{\ast }}{a\alpha \gamma }(\frac{A^{\ast }(a{A^{\ast }}^2-b)}{A_p^{\ast }}+2b)+2d\beta A^{\ast }](\alpha L^{\ast }+d{A}^{\ast })\\ & -\frac{\beta d{A^{\ast }}^3}{A_p^{\ast }{A}_c^{\ast }}(\alpha L^{\ast }-d{A}_p^{\ast })-\frac{2\beta d^2{A^{\ast }}^3}{A_c^{\ast }}.\end{array}$ This indicates the following two cases:

1. If Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has a unique interior equilibrium $E_2=(L_2,A_2,{A_p}_2)=(\frac{d}{\beta }{A}_2^2,A_2,\frac{\alpha \beta A_2+\alpha A_2^2}{\beta +\alpha A_2})$, then under the conditions $a>\frac{bd(1-\alpha )}{{\alpha }^2\gamma }$, $a_1<a<a_2$ and $\begin{array}{rl}& max\{0,\frac{\beta }{\beta -\alpha A_2}\}\\ & <\frac{A_2}{{A_c}_2}<min\{\frac{\alpha L_2+d{A}_2}{\alpha L_2+d{A_p}_2}[\frac{d(a{A_2}^2-b)}{a\alpha \gamma }+\frac{2{A_p}_2}{A_2}(\frac{bd}{a\alpha \gamma }+1)],M\},\end{array}$ where (17) $\begin{array}{rl}{a}_1=& \frac{b({\alpha }^2\gamma +4d{\beta }^2(1-\alpha )-\alpha \sqrt{\gamma ({\alpha }^2\gamma +8d{\beta }^2(1-\alpha ))})}{2{\alpha }^2{\beta }^2\gamma },\\ a_2=& \frac{b({\alpha }^2\gamma +4d{\beta }^2(1-\alpha )+\alpha \sqrt{\gamma ({\alpha }^2\gamma +8d{\beta }^2(1-\alpha ))})}{2{\alpha }^2{\beta }^2\gamma },\\ M=& \frac{1}{{\beta }^2}(2\alpha d{A_2}^2+(\alpha +2d)\beta A_2+3{\beta }^2)+\frac{d}{a\alpha \gamma }(\frac{A_2}{{A_p}_2}(a{A_2}^2-b)+2b),\end{array}$(17) we get $c_2>0,c_3>0$. And we can verify that $c_1{c}_2-c_3>0$. Thus, we can conclude that the interior equilibrium $E_2$ is locally stable by applying the Routh-Hurwitz criterion.

2. If Model (4(4) $\begin{array}{rl}{L}^{\prime }& =\gamma S_{L_{max}}\frac{a{A}^2}{b+a{A}^2}-\beta L,\\ A^{\prime }& =\beta L-d{A}^2,\\ A_p^{\prime }& =\beta L\cdot max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}+(max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_p}{A_c}}\}{A}_c-max\{0,\frac{\mathrm{\partial }{S}_L}{\mathrm{\partial }\frac{A_c}{A_p}}\}{A}_p)L-dA{A}_p.\end{array}$(4) ) has two interior equilibria $E_i=(L_i,A_i,{A_p}_i)=(\frac{d}{\beta }{A}_i^2,A_i,\frac{\alpha \beta A_i+\alpha A_i^2}{\beta +\alpha A_i}),i=1,2$ where $E_1<E_2$, then under the conditions $a^{\ast }<a<\frac{bd(1-\alpha )}{{\alpha }^2\gamma }$, $a_1<a<a_2$, and $\begin{array}{rl}& \frac{A_1}{{A_c}_1}>\frac{1}{{\beta }^2}(2\alpha d{A_1}^2+(\alpha +2d)\beta A_1+3{\beta }^2)-\frac{d}{a\alpha \gamma }(\frac{A_1}{{A_p}_1}(b-a{A_1}^2)-2b),\\ & \frac{A_2}{{A_c}_2}<\frac{1}{{\beta }^2}(2\alpha d{A_2}^2+(\alpha +2d)\beta A_2+3{\beta }^2)+\frac{d}{a\alpha \gamma }(\frac{A_2}{{A_p}_2}(a{A_2}^2-b)+2b),\end{array}$ we obtain $A_1<\sqrt{\frac{b}{a}}<A_2$ and $c_2(E_1)<0$ but $c_2(E_2)>0$. We also can verify that $c_1(E_2)c_2(E_2)-c_3(E_2)>0$. Therefore, the interior equilibrium $E_2$ is locally asymptotically stable while $E_1$ is unstable.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research of F.R. is partially supported by the National Science Foundation of China [grant numbers 11601226, 11426132 and 71871115], Qing Lan Project of Jiangsu Province, the Natural Science Foundation of Jiangsu Province of China [grant number BK20140927], and the research funds from Nanjing Tech University and Jiangsu Government Scholarship for Overseas Studies. The work of Y.K. is also partially supported by NSF-DMS [grant numbers 1313312 and 1716802]; NSF-IOS/DMS [grant number 1558127], DARPA – SBIR 2016.2 SB162-005 Phase II, and UHC Scholar Award 220020472.