Nutritional Regulation Influencing Colony Dynamics and Task Allocations in Social Insect Colonies

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 labor 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.


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 behavioral connections between its members [1,2]. They exhibit a decentralized system with a sophisticated division of labor resulting from interactions among members of the colony and the environment [2][3][4]. In addition to the reproductive division of labor between the queen and the workers, workers also have a division of labor between foragers which leave the nest to search for food and non-foragers that carry out tasks within the nest [5]. In social insect societies, foraging responsibilities are assigned to a subset of adult colony members [6]. Internal and external factors happening at both the individual and colony level shape the foragers' decision to bring back a certain type of food [6].
Currently, there are few studies that have focused on the outcome of nutrient regulation in social insects at the colony level [6] (but see [7][8][9]). 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 behavior 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 indivdiual forager is based not only on the colony's current nutritional status, but also on the worker's physical caste, age, and prior experience [10,11]. The nutritional needs of the colony are shaped by the differing needs of larvae and workers in the colony [6,11,12]. For instance, the growth of larvae relies heavily on protein, while worker ants require primarily carbohydrates as a source of energy [12][13][14][15][16][17]. Many studies have shown that the ratio of protein to carbohydrates in the diet of a range of insect species is crucial for performance [6,12,[18][19][20], though, in general, carbohydrates are often more attractive to foragers than protein [7,11]. However, the protein required for growth may be in greater demand when a queen is laying eggs [11].
In order for social insect foragers to compensate for potential nutrient restrictions in the food available to the colony [11,[21][22][23], foragers adjust their collection in favor of food sources containing limiting nutrients [7,11,23]. This guarantees that the colony meets its longer term objectives and thus promotes colony growth and reproduction [11]. According to Dussutour et al. [7], within a colony, workers recruit nestmates for food collection at different rates depending upon food type [16,24], food concentration, and hunger level [7,25]. At the individual level, when workers are starved recruitment will be stronger to carbohydrate-rich food sources than to sources high in protein [7]. 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 [7,24,26].
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 [6-8, 17, 27, 28]. However, currently there are no mathematical models to our knowledge that have attempted to study these mecha-nisms 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 analyze 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 labor implicitly. Also, by considering the basic mechanisms affecting colony growth such as cooperative effort for reproductive division of labor, 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.

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 : The brood population L increases with the average egg-laying rate of the queen(s) given by γ, which is discounted by two factors: (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. [29,30], where the cooperative efforts that lead to the eggs' survival is measured by a Holling type-III function aA 2 b+aA 2 , where b is a half-saturation constant and a is the portion of the division of labor invested towards the successful development of the larvae. (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 Ap Ac .
Examples of S L could be S L Ap Ac = −α 1 Ap Ac − θ m + α 2 (θ c − θ m ) with α i ∈ (0, 1) i = 1, 2 as a scaling factor of nutrient collection, θ m representing the optimal nutrient ratio, and θ c representing the maximal nutrient ratio that brood can survive (see Figures 1(a)), or general functions such as the normal biological performance curve (see Figure 1(b)). Notice that S L ≤ 1 can be negative, thus we define S Lmax = max{0, S L } such that S Lmax ∈ [0, 1] is a survival probability.
The brood population decreases by a maturation rate βL, which describes the rate at which brood matures into the adult class A. Thus, we have following equation:

When
Ap Ac is less than θ m , the brood survival rate increases, and decreases when Ap Ac is greater than θ m . This phenomenon has been supported by the work of [6,8,20,31], 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 In this study, we assume that when the nutrition level hits or exceeds the critical value θ c , i.e., Ap Ac ≥ θ c , the nutrient becomes toxic such that no brood can survive. Lastly, Figure 1   In general, we expect the protein to carbohydrate ratio Ap Ac of the colony to fall in a certain range in order for the colony to survive and grow, say, Ap Ac ∈ [θ 0 , θ m ], where θ 0 is the minimum nutrient ratio necessary for brood survival. This is supported by [6,8]. Thus it is reasonable to assume that S Lmax Ap Ac has the following simple form: subject to α i ∈ (0, 1), i = 1, 2 and α 1 (θ m − θ 0 ) + α 2 (θ m − θ c ) = 0 and 0 < α 1 (θ m − θ 0 ) ≤ 1.
In particular, we have Then we have In the symmetric case, i.e., α 1 = α 2 = α, then we have Ac < θ c . In addition, we have the following: The special case of the symmetric scenario above is S L (0) = 0, i.e., θ c = 2θ m . In this case, we have S L Ap Ac = −α Ap Ac − θ m + αθ m with θ m ∈ (0, 1 α ) in our proposed model (4). In the following section, we will provide mathematical analysis of the general case of S Lmax Ap Ac shown in (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 dA 2 . The density-dependent mortality rate follows the approach of [29], thus we have following equation: average mortality rate

The protein forager population A p (t): The ratio
Ap Ac measures the nutritional status of the colony. Assume that the brood can survive in a range of nutrient ratio, i.e.

Ap
Ac ∈ (θ 0 , θ c ), and any ratio greater than θ c can be toxic to brood. In addition, there is an optimal nutritional ratio Ap Ac , denoted by θ 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 Ap Ac when Ap Ac ∈ (θ 0 , θ m ), and passing this optimal ratio θ m , the brood survival rate decreases with Ap Ac . The survival rate of brood is zero when • The portion of the successful brood developed into adults which enter into the protein forager pool can be modeled by the term: represents the nutritional requirements of the colony measured by the ratio Ap Ac 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 [7,11,16,[23][24][25]. 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 Ac is less than the optimal nutrient ratio θ m (where the maximum brood survival rate occurs), then we expect ∂S L ∂ Ap Ac > 0 and ∂S L ∂ Ac Ap < 0, thus this indicates that carbohydrate foragers will switch tasks to forage for protein, i.e. max 0, ∂S L ∂ Ac In a similar fashion, if the ratio θ m ≤ Ap Ac ≤ θ c , protein foragers will switch to forage for carbohydrates, that is, max 0, • The total forager population A = A p + A c decreases with a density-dependent death rate dA 2 = dA(A p + A c ), then the protein forager population decreases with the densitydependent mortality rate dAA p .
Considering the factors above, we derive the population dynamics of the protein forager as follows: portion of matured adults entering Ap Based on the ecological assumptions above, the population dynamics of a social insect colony with nutrient regulating foraging activities is described as follows: The biological meaning of the parameters and the related values are listed in Table 1.

Mathematical analysis
The state space of the proposed ecological model (4) is R 3 + . All parameters a, b, d, α, β, γ, θ 0 , θ m , θ c are assumed to be strictly positive based on their biological meaning. We focus on the proposed function S L Ap Ac shown in Eq. (1) and Figure 1(a). The related mathematical results should be easily adopted to the symmetric case (3). Under such conditions, we first show that Model (4) is biologically well-defined, i.e., it is positively invariant and bounded in R 3 + in the following lemma: The extinction equilibrium E 0 = (0, 0, 0) of Model (4) always exists. The local stability of the trivial equilibrium E 0 cannot be analyzed directly for our model (4). However, from the first two equations of Model (4), we have Note that Model (4) is positively invariant and bounded from Lemma 3.1, thus we can conclude that if aα 1 γ(θm−θ 0 ) b < d, then the inequality above implies that lim sup t→∞ (L + A) converges to a nonnegative constant. In addition, we have L | A=0, (4) converges to the extinction equilibrium E 0 in R 3 + . Thus, we have the following proposition: Remarks: Note that the inequality . 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 θ 0 = 2θ m − θ c , then the inequality becomes aαγ(θc−θm) b < d with α 1 = α 2 = α.
Assume that E * = (L * , A * , A * p ) is an interior equilibrium of Model (4)  Thus, we have To solve for (L * , A * , A * p ), we set L = A = A p = 0, which implies the following equations Therefore, A * of an interior equilibrium (L * , A * , A * p ) satisfies the following equation: Recall that α 1 ∈ (0, 1). Depending on the exact values of a, b, d, α 1 , β, γ, θ 0 , θ m , θ c , the equation (5) can have zero, one, or two positive roots. Let be the possible positive roots of equation (5). Let us denoteâ * as follows: which is an increasing function of θ 0 .
In the symmetric case, we have θ 0 = 2θ m − θ c , thenâ * shown in (6) can be rewritten as Also note that Then the following theorem provide conditions for existence of equilibrium solutions of Model (4): , then there is only one trivial equilibrium E 0 = (0, 0, 0) and no other positive equilibrium.

If a >
bd(1−α 1 ) α 1 γ(α 1 −(1−α 1 )θ 0 ) and θ 0 < α 1 1−α 1 , then Model (4) has only one positive equilibrium E 2 1−α 1 , then Model (4) has two positive equilibria in the following form in addition to E 0 : Remarks: The detailed proof of Theorem 3.1 is shown in the last section. The number of equilibria of Model (4) is determined by the positive root(s) of equation (5). Theorem 3.1 implies that the value of the division of labor invested on larvae a and the minimal protein to carbohydrate ratio θ 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) has simple dynamics: no limit cycle and only equilibrium dynamics. At the stable equilibrium, the ratio describing the nutritional level of the colony is Equation (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 * depends on θ 0 , so does as well.
Now we discuss stability of the interior equilibrium for Model (4). Let E * = (L * , A * , A * p ) be an arbitrary positive interior equilibrium of Model (4). The Jacobian matrix associated to Model (4) at equilibrium is: Then the characteristic equation of J| E * is where C 1 = β + αL * + 3dA * > 0, The stability of the steady state E * = (L * , A * , A * p ) can be determined by the distribution of the roots of Eq. (10). That is, if all the roots of Eq. (10) have negative real parts, then E * is locally asymptotically stable; if at least one root of Eq. (10) has positive real parts, then E * is unstable; if any root has zero real part and other roots all have negative real parts, then the stability of E * cannot be determined by the linearized system directly.
The following theorem provides a global result on dynamics of the proposed model (4) regarding when a colony will collapse. Theorem 3.2 (Extinction of species). If 0 < a <â * and θ 0 < , then Model (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 labor invested on larvae a and the nutrient θ 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.
Biological Implications: The results in Lemma 3.1, Theorems 3.1 and 3.3, imply that the division of labor invested on larvae a decreases past the critical pointâ * = shown in (6) and the first dotted line in Figure 2. Model (4) exhibits a backward bifurcation shown in Figure 2 where b = 0.1, d = 0.1, α 1 = 0.3, β = 0.7, γ = 0.9. In Figure 2(a), we set θ 0 = 0.1 and in Figure 2(b), we set θ 0 = 0.2. Based on the expression of the critical valueâ * shown in (6),â * is an increasing function of θ 0 , which is reflected in the difference between Figure 2(a) and Figure 2(b). The value of θ 0 measures the minimum ratio of protein to carbohydrates that can allow the survival of larvae. Simulations shown in Figure 2(a), 2(b) and 3 suggest that the larger value of θ 0 , the more likely the colony can survive with a larger population of workers A and thus the higher nutrient ratio Ap Ac . In summary, our theoretical work combined with the related simulations suggest that the division of labor invested on larvae a and the minimal nutrient ratio θ 0 can affect colony survival, the distribution of the brood and the total forager population affect the protein forager population. For instance, the larger θ 0 , the more division of labor 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 Ap Ac versus the division of labor invested on larvae a with different values of θ 0 , other parameters values are taken as those in Figure 2. 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â * 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 color indicates E2 which is always stable; the green color indicates E1 which is always unstable; and the red color is the extinction equilibrium E0.
Remarks: Theoretical results and numerical simulations (see Section 4) confirm that the special symmetric case of Model (4), i.e., θ 0 = 0 and θ c = 2θ m , undergoes a backward bifurcation as a decreases past the critical value a * defined in (12). Conditions shown in Theorem 3.4 suggest the importance of the protein to carbohydrate ratio, i.e., A 2 Ac 2 = 1 + Ap 2 Ac 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) undergoes a backward bifurcation as a decreases pastâ * (or a * in the case of θ 0 = 0 and θ c = 2θ m ). More specifically, it exhibits the following global dynamics: 1. If 0 < a <â * and θ 0 < , then the colony collapses due to the lack of efforts of division of labor invested on larvae and the minimum nutrient requirement θ 0 being too low.
, the survival of the colony depends on its initial population size.
Our theoretical results suggest that the survival rate of larva to worker S Lmax Ap Ac plays critical roles in determining colony population dynamics. We assume that S Lmax Ap Ac takes the form of (1) based on relevant biological studies. Our analysis implies that the values of parameters α 1 and θ 0 in S Lmax Ap Ac have pronounced impacts on dynamical outcomes. In the next section, we use bifurcation diagrams to explore detailed impacts.

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 Lmax Ap Ac , the dynamics of Model (4) depends on the division of labor a, egg laying rate γ, the scaling factor on the brood survival rate due to the nutritional status α 1 , the minimal nutrition ratio θ 0 , the maturation rate β, and the natural mortality d. To explore the effects of a and θ 0 , we perform bifurcation diagrams in Figure 2   We can see that the larger value of a can lead to the larger population A (see Figure 2) and the larger nutrient ratio Ap Ac (see Figure 3). Figure 2 and 3 also show the effects of the minimal nutrient requirement for brood survival θ 0 : The larger value of θ 0 , (1) the larger critical thresholdâ * ; (2) the smaller population A; and (3) the smaller nutrient status, i.e., the smaller value of Ap Ac . Next, we perform bifurcation diagrams of Model (4) regarding how the minimum nutritional requirement θ 0 and the scaling factor of survival probability of brood α 1 affect population dynamics of the colony in Figure 4. Figure 4(a) shows that Model (4) exhibits reversed backward bifurcation on θ 0 . Figure 4(b) suggests that the larger value of α 1 , the larger population of worker A and the better probability of colony survival.
Symmetrical case θ c ≤ 2θ m (i.e., θ 0 = 2θ m − θ c > 0): Figure 5(b) provides an example of bifurcation diagram on division of labor invested on larvae (a) of Model (4) by taking same parameter values in Figure 5(a) except that θ c = 7.8. Figure 5(b) shows that Model (4) undergoes a bifurcation as the portion of the division of labor invested on larvae a decreasing pastã * =

Conclusion
Variation in nutrient consumption among individuals is considered a conserved mechanism regulating castes and division of labor 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 [32]. 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 favor of food sources containing limiting nutrients to maintain colony growth and reproduction [7,11,23].
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 Lmax Ap Ac . Our formulation of S Lmax Ap Ac is based on biological studies (see [7,11]) and embeds with an adaptive modeling approach adopted from [29,33,34]. 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) can be summarized as follows:

The nutrition status measured by
Ap Ac 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 Ap Ac , 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: (a) The division of labor 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) exhibits a backward bifurcation (shown in Figure 2) as a decreasing past its critical point which is an increasing function of the maturation rate β, the minimal nutrient ratio θ 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.  Figure 4(b) suggests that the larger value of α 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 behavior of leaf-cutter ants during its ontology.