Effect of interaction network structure in a response threshold model

Task allocation is a characteristic feature of social insects. This phenomenon is self-organized by workers in a colony without central instructions. Workers perform the necessary tasks while communicating and determining their colony’s local situations or outside environments. They can also allocate the workload of the task to finish quickly. To explain the self-organized phenomenon, we develop a theoretical model that includes the local interaction mechanism among the workers. We extend the fixed response threshold model by introducing the interaction network among workers. The new model is called the interaction network response threshold model. This model analyzes the effects of the interaction network structure on the workload allocation mechanism using the Gini coefficient. We find that the network structure affects the information diffusion process, and the network features affect workload allocation. Our results indicate that communication among the workers in a colony is an essential system for workload allocation.


Introduction
Task allocation is a characteristic feature of social insects, such as ants and honeybees [6]. A colony of ants comprises hundreds of individuals, and the so-called workers have to perform numerous tasks to survive. These tasks are variational and massive, including feeding, keeping nests clean, and caring for their larvae. In addition, the necessary amount of workload for the tasks varies from time to time. Therefore, effective task allocation mechanisms are essential to keep an ant colony healthy. Despite the significant number of tasks requiring an effective task and workload allocation mechanism, the workers assess the situation based on local information and perform tasks with an appropriate amount of workload through self-decision without leaders in a colony. To explain this mechanism, Robinson proposed the response threshold hypothesis [12]. Mathematical models for the hypothesis [2,15] explain the mechanism of task allocation. Particularly, in the fixed response threshold (FRT) model proposed by Bonabeau et al. [2], workers stochastically decide whether to perform their tasks, and all workers can precisely access the global stress. The FRT model succeeded in describing how different tasks were allocated based on the difference in the workers' thresholds for the specific tasks by quantifying the workers' activities for each task. In addition, the model could explain the differentiation between lazy and diligent workers for a task [7]. The FRT model has the ability to explain the workload distributing mechanism. However, these assumptions differ from the natural state of the ant colony and have not been experimentally confirmed. Furthermore, the local interaction network in ant colonies is not considered in the FRT model, although numerous studies have reported that the network is essential to organizing the colony through task allocation [8-11, 13, 16]. Greenwald et al. [4] presented the method visualizing the interaction of trophallactic as a network in a colony, and Greenwald et al. [5] quantified the food distribution in a colony with entropy.
In this study, we introduced the interaction network structure into the FRT model; the new model is called the interaction network response threshold (INRT) model. In the INRT model, individual workers share local stress information through interaction based on the interaction network. We adopted three network structures: two-dimensional lattice, scale-free, and small-world network structures. We analyzed the effects of the different network structures for the interaction network on the workload allocation mechanism using the Gini coefficient, which can quantify the inequality in economic wealth distribution [1]. Consequently, the local interactions in the two-dimensional lattice disturb the development of labor division during workload allocation. In contrast, the local interactions in a scale-free and small world enable the development of labor division. Finally, we report the relationship between the detailed network properties and the workload allocation phenomenon. All simulations were conducted using Python 3.9 with the NetworkX 2.6, NumPy 1.16, and Matplotlib 2.2 packages.

Model and analysis
We introduce the local interaction among workers into the FRT model. To explain our model, we review the FRT model and explain the details of our new model, the INRT model. Furthermore, we introduce the Gini coefficient and represent the analysis method of the dynamics of the FRT and INRT models.

Fixed response threshold model
The FRT model proposed by Bonabeau et al. [2] is a mathematical realization of the response threshold hypothesis proposed by Robinson [12]. The FRT model has two essential hypotheses: all workers share the stress for the entire colony, and workers engage in tasks without communicating with their coworkers, as shown in Fig. 1a. The FRT model consists of N individuals, and X i (t) is the state of each individual, where i indicates the individual ID. X i (t) can be either an active state ( X i (t) = 1 ), wherein the individual takes action to reduce stress S(t), or a resting state ( X i (t) = 0 ), wherein the individual does nothing. The state X i (t) is probabilistically chosen by the function of stress S(t) and individual response threshold i , where the stress S(t) on the entire system evolves in time according to the following equation: where is the stress increase rate, is the stress decrease rate per worker, and N active (t) = ∑ N i X i (t) . The probability of ants changing their state is defined by the following two equations: where p c is the constant value corresponding to active ants that quit working. The FRT model comprises the aforementioned three equations; furthermore, the system is defined by the parameters N, , , and p c and the distribution of i . (1)

Local stress diffusion
We introduce interactions among workers to share stress information. In contrast to the FRT model, limited workers can assess stress S(t) in the INRT model. Each worker has local stress information i (t) instead of S(t), and i (t) is determined based on the following equation: where k i represents the edge degree of the interaction network G , and M ij is the component of the adjacent matrix M of G , so that M ij = 1 between workers i and j if interaction exists, and M ij = 0 otherwise. We define the interaction network G as an undirected graph, because we do not consider the direction of information transfer. The interaction network G is important, because it directly determines the interaction of workers in simulation. The detail of how to construct the networks will be explained in the next subsection. In the INRT model, information of stress S(t) diffuses into the entire colony based on Eqs. (2) and (3) through the interaction network G . The probability of changing the state of workers in the INRT model follows the same function as described in Eqs. (2) and (3)

, but S(t) is substituted by i (t):
In this study, we define a k-th worker that can assess stress S(t) for simplicity, as shown in Fig. 1b. The worker exceptionally determined her behavior based on S(t); thus, k (t) = S(t) . The worker that can assess the stress was randomly chosen in every simulation without considering its degree.

Interaction network structures
The interaction network structures define the dynamical aspects of local stress information. We consider three types of interaction networks: a two-dimensional lattice network, a scale-free network, and a small-world network. The three interaction networks were constructed using functions in NetworkX 2.6 explained in Table 1 and converted to the corresponding adjacent matrixes. The two-dimensional latticelike network was generated, so that workers have edges as if they were allocated on a two-dimensional finite lattice, as shown in Fig. 2a. The workers inside the lattice, the workers on the edge, and the workers on the corners have four connections, three connections, and two connections, respectively. The average degree of the two-dimensional latticelike network < k > equals to 4(1 − 1∕ √ N) , thus < k >= 3.6 when N = 100 . The scale-free network structure was generated as a random graph based on Barabási-Albert (BA) model (Fig. 2b). The scale-free network generated by BA model has the degree distribution which satisfies a powerlaw distribution, P(k) ∼ k − , where k means the degree of  The third column is the average degree for each NW structure. N is the number of nodes that represents the number of workers. m is the number of edges to attach from a new node to existing nodes. k ws defines that each node is joined with its k ws nearest neighbors in a ring topology. p is the probability of rewiring each edge NW structure Generating function Scale-free barabasi_albert_graph(N, m) 3.92 Small-world watts_strogatz_graph(N, k ws , p) 2.0 each node and ≈ 3 in BA model. In the case N = 100 , < k >≈ 3.92 . The small-world network structure was also generated as a random graph based on Watts-Strogatz (WS) model (Fig. 2c). In this study, we defined k WS = 2 and p = 1 , and the average degree became that < k >= 2.0 . We compare the extended model with these network structures and the conventional FRT model using Gini coefficients. In the next section, we describe the comparison method in detail.

Workload distribution
We define the workload of each worker, A i , using the following equation: The workload distribution is calculated by the rank distribution of A i . We rearranged the index i of A i such that A i < A j for i < j . In this study, A i < A j for i < j is satisfied for all results. Subsequently, the relative workload a i , which is the relative ratio of each worker's workload to the total workload in the colony, is defined by the following equation: We calculate the relative workload to compare the results with different total workloads.

Gini coefficient
The two models explain the dynamics of task allocation in social insects and generate a workload distribution in which a few active workers did the task. Therefore, we calculate the Gini coefficient, G, which demonstrates the degree of inequality in income distribution, to quantify the inequality of the relative workload distribution [1]. The Gini coefficient, G, is defined by the below equation All workers become an active state with equal frequency when G = 0 . On the other hand, a few workers contribute to reducing the stress, and the others were inactive states when G = 1 . Figure 4 shows the phase diagram of G for the simulation results of the FRT model. To explain the task allocation phenomenon on the workload distribution, the parameter condition, > , must be satisfied. That condition makes the value of G larger. If , the stress increase rate, is too strong ( < ), all workers have to work as many as possible, and thus, G becomes a smaller value.

Numerical results
We simulated the colony activities for the FRT and INRT models using multiple parameters. Subsequently, we calculated the relative workload distributions and the Gini coefficient of the distributions to evaluate their inequalities, described as follows. First, we calculated the Gini coefficient G of the relative workload distribution with different values of , , and p c in the FRT model. Next, we calculated the Gini coefficients with the same parameters mentioned above in the INRT model with different interaction network structures.
The following parameters were common to all simulations: the number of workers N is 100, simulation time T is 2000, the threshold i for each worker follows a normal distribution N(1, 2), and a number of simulation ensembles are 100. We confirmed that the simulation time is enough by observing that the time series of A(t) = ∑ k A k (t) and S(t) become stable. The analysis of these time series is our future work. As for i , we take the absolute value of i , if i was a negative value. Initial conditions for the variables were defined, so that i (0) = 0 , X i (0) = 0 for all i, and S(0) = 0 for all simulations.

FRT model
In the FRT model, the three parameters , , and p c in Eqs.
(1)-(3) determined the workers' activity and generated the relative workload distribution with a specific shape. Figure 3a, b shows a i when p c = 0.1 and = 0.1 and = 0.5 , respectively. As increased, the relative workload distribution became unequal. Figure 3a, b shows that where < , each worker reduces a relatively small amount of stress; therefore, a large proportion of workers has to be active. Where > , the small proportion of workers with low reaction thresholds took most of the workloads, because individual workers reduce a large amount of stress. Therefore, the workloads of workers with small reaction thresholds increased, whereas the workloads of other workers with high reaction thresholds did not increase, which resulted in inequality in the workload distribution. A division of labor was established because of the dynamics by and . The phase diagram of the Gini coefficient G FRT for the combination of and with p c =0.1 shows how the division of labor is established for each parameter.

INRT model
Next, we compared the simulation results of the three network structures shown in Fig. 2 incorporated into the interaction network with the results of the FRT model obtained 1 3 in the previous section. The Gini coefficients for the INRT model in which the two-dimensional lattice network, scalefree network, and small-world network structure were introduced are denoted as G lattice , G scale , and G small , respectively. The phase diagram of G FRT for the FRT model (Fig. 5c) and the phase diagrams for the Gini coefficients G lattice , G scale , and G small are qualitatively similar, as shown in Fig. 5.
To quantitatively compare the results of the FRT model with the ones of the INRT models, the phase diagrams of G lattice − G FRT , G scale − G FRT , and G small − G FRT for the same values of stress increase rate and task throughput are shown (Fig. 6).
The phase diagrams for the three differences, G lattice − G FRT , G scale − G FRT , and G small − G FRT , differ. This indicates that the structure of the relative workload distribution is affected by the interaction network.
First, we discussed the trend common to all models that introduced the network structure. The Gini coefficients of the INRT models were smaller than G FRT in the parameter region > , where the colony established the division of labor. The reduction in the Gini coefficient indicates that the interaction network weakens the unequal activity distribution. In the FRT model with a parameter where the > , the workers with low reaction thresholds immediately reacted to the increase in stress and became active in decreasing stress, which made it difficult for the other workers with high reaction thresholds to become active. The workers in the FRT model could assess the stress state; therefore, the colony could quickly respond to an increase in stress. However, the stress does not propagate rapidly in the entire colony, unlike the case in the FRT model where the workers share stress information through the interaction network following the diffusion effect in Eq. (4). Because of the slow propagation of the stress state in the colony, workers with a low response threshold could not detect the stress increment, even though the stress was high enough to make the workers active. In contrast, the high-stress state in the INRT model relative to the one in the FRT model made the workers with a highresponse threshold active more frequently. Consequently, the above mechanism could weaken the skew relative workload distribution. In the following sections, we discuss the changes in the Gini coefficients for each network structure in the INRT model from the ones for the FRT model separately.

Two-dimensional lattice network
The Gini coefficients for the INRT model with the twodimensional lattice interaction network, G lattice , decreased by a maximum of approximately 0.1 for p c = 0.1 (Fig. 6a) in the region where > . The reduction in G lattice indicates that the model had a weaker division of labor structure than that in the FRT model. This tendency of reduction in the Gini coefficient was similar when p c = 0.5, 0.7 (Fig. 6b, c), but the maximum decreases were less than 0.1. In addition, when the rest probability p 0 was high, the interaction network affected the system's dynamics to a low degree.
The slight reduction was because i (t) takes a small value in the INRT model. i (t) is smaller than S(t) because of the diffusive effect by Eq. (5). If i (t) is small in Eq. (5), the rest parameter i primarily determines the probability P(X i ∶ 0 → 1) . Therefore, the response threshold distribution was important to determine the relative workload distribution. The response thresholds of the workers were random variables with a Gaussian distribution; therefore, the relative workload distribution was not unequal like in the FRT models.

Scale-free network
In the case wherein the interaction network had a scale-free network structure, the reduction in the Gini coefficient was similar to the results with the two-dimensional lattice network structure for > when p c = 0.1 (Fig. 6d). In contrast, the Gini coefficients increased when ≤ , which indicates a strong tendency for inequality of the relative workload distribution under stressful environments. However, the values of the Gini coefficients were too small to consider that these changes were relative to the FRT model and that a division of labor was established (Fig. 5d-f). However, the values of G scalefree were greater than those of G FRT and were not large enough to consider that the workers generated division of labor.

Small-world network
In the case wherein the interaction network had the smallworld network structure, the relative workload distribution of the workers was more equal than that in the FRT model. In particular, when p c = 0.1 (Fig. 6g), the Gini coefficient decreased by approximately 0.2, thereby indicating a fairly strong homogenizing effect. In contrast, in other regions, the Gini coefficient increased by approximately 0.1 (Fig. 6g-i).
Compared to the scale-free network structure results, the relatively large increase in the stressful environment is a remarkable feature in the results of the small-world network structure.

The results' comparison
Finally, we discuss the similarity in the scale-free network and small-world network structure results. The G scalefree and G small have similar phase diagrams for and . One of the possible reasons is that the scale-free network has the features of a small-world network and vice versa. Essentially, the similarities in the results of G scalefree and G small reflect the intrinsic commonalities among the two network structures.

Summary and discussion
We analyzed the reaction threshold hypothesis, which explains the division of labor in communities of social insects, such as ants and bees, and its mathematical model, the FRT model. In addition, we highlighted the discrepancy between the FRT model and the latest observations, proposed an extended model, the INRT model, and discussed their possibilities. By introducing a specific network structure in the interaction of workers, we found that the relative workload distribution has a specific structure even if workers behaved based on the same dynamics. In a two-dimensional lattice network, where information is exchanged only with neighboring workers, workers are less likely to establish a division of labor than in the FRT model, where they can measure the stress state directly. Information exchange by interaction networks with a small-world or scale-free network structure, rather than a two-dimensional lattice network structure, showed that workers might be able to propagate their stress state to all workers more efficiently. This result suggests that the workers' intrinsic reaction thresholds and stress-induced dynamics do not determine the worker's activities. Nevertheless, the information exchange among workers has a vital role in establishing the division of labor. In this study, we only considered a limited number In addition, as shown in previous studies [8,9], network structure has been experimentally confirmed, and its relationship with the change in worker roles has been clarified. Swain et al. [14] analyzed the experimental data by Mersch et al. [9] in more detail, and showed that the network structure was strongly related to the dynamical development of task allocation in a colony. Especially, they reported that workers that changed their role in a colony have a higher degree as a network node. In the case of the existing data by [9], the number of nodes, the workers, was from 40 to 160, and the workers interacted with almost 70% of their nest mates; the average degrees were very high. The network features in our theoretical research were far from the ones in the observed data. It is important to simulate the INRT model with parameters that are close to observed data, but we have to pay attention to the difference in the time scales in our theoretical model, and the observed data. The existing data aggregated for a day unit included many kinds of behaviors in a colony. In contrast, we focused on the workload distribution on one task; thus, effective time scale of the FRT model has to be hours. Therefore, we have to analyze the NW features on a shorter time scale experimentally to validate the theoretical model we proposed. To specify the function of the interaction network in the task allocation mechanism, we have to observe and analyze specific dynamical behaviors like foraging or nursing.
In the future, we aim to analyze the interaction between workers through image analysis using fiducial markers (ex. ARTag [3]). The network structure parameters obtained from the experiments will be introduced into the INRT model to verify the model's validity more precisely. In addition, we observed the activity time series with temporal group structures obtained by our research group in an experiment on ants using RFID tags [17]. Because no mathematical model has explained this activity time series, our task is to clarify the relationship between the network structure in the colony and the activity time series using the INRT model.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.