Animals live in groups to defend against predation and to obtain food. However, for some animals—especially ones that spend long periods of time feeding—there are costs if a group chooses to move on before their nutritional needs are satisfied. If the conflict between feeding and keeping up with a group becomes too large, it may be advantageous for some groups of animals to split into subgroups with similar nutritional needs. We model the costs and benefits of splitting in a herd of cows using a cost function that quantifies individual variation in hunger, desire to lie down, and predation risk. We model the costs associated with hunger and lying desire as the standard deviations of individuals within a group, and we model predation risk as an inverse exponential function of the group size. We minimize the cost function over all plausible groups that can arise from a given herd and study the dynamics of group splitting. We examine how the cow dynamics and cost function depend on the parameters in the model and consider two biologically-motivated examples: (1) group switching and group fission in a herd of relatively homogeneous cows, and (2) a herd with an equal number of adult males (larger animals) and adult females (smaller animals).

Although animals gain many advantages—such as protection from predators—from living in groups, they also incur considerable costs. For grazing animals, such as cows and antelopes, these costs include balancing their own nutritional needs to stay in one place to feed with the need to keep up with a group and stop grazing when the rest of the herd moves on. If the nutritional needs of different individuals are sufficiently disparate, this can lead to the splitting of a group so that those with similar needs to graze, lie, and ruminate remain together. If a group of animals becomes too small, however, this can increase the risk of predation, as small groups are less able to defend themselves against predators than large groups. In this paper, we describe a cost function (CF) that balances predation risk (based on group size) with different individual needs for feeding and lying down to infer the sizes at which group splitting occurs. We model variation in hunger and lying desire using the standard deviation of individuals within a group, and we model predation risk as an inverse exponential function of the group size. In a series of examples, we optimize the CF for each individual in a group of animals and examine when groups of cows split into smaller groups.

## I. INTRODUCTION

Animals gain many advantages from grouping and synchronizing their behavior—including greater vigilance, coordinated defense against predators, and increased ability to find and defend food sources.^{1,2} However, living in large groups also carries disadvantages, such as increased risk of disease and parasitism,^{3,4} having food stolen,^{5} and interference with movement.^{6} A “perfect” synchronization requires animals to change their activities at a communal time rather than at individual ideal times, and this can be costly for individuals.

The balance between synchrony and risk of predation is complex,^{7,8} and one possible approach for examining such a balance is with a cost function (CF) with components from synchrony and risk. When a group of animals becomes very large, the cost incurred through synchrony tends to exceed that incurred through risk, as a significant number of individuals change their desired activities (like eating or lying) to conform to communal decisions. Because of the balance, a CF with components from synchrony and risk of predation should have at least one optimum point, and one should expect animal groups to split if they are too large. However, an optimal group size is not necessarily a “stable” group size. Supposing animals join a group one by one, a stable group size is a size at which there is no further fission of groups or switching of animals between groups.^{9} Even when a group is already at its optimum for existing individuals, extra individuals can still benefit from joining the group. At some point, however, the group can become sufficiently large that it splits into two groups, as this benefits its members more than the overloaded single group.^{9} A stable group size is therefore likely to be consistently larger than an optimum group size.^{9,10}

Sometimes grouping can be even more complicated, as individuals within a group differ in many ways that relate to their fitnesses. For example, males and females in a herd differ in their nutritional needs. However, although they can benefit from staying with a mixed-sex group, some individuals may have to interrupt valuable feeding or lying time to keep up with a herd when it moves.^{11,12} It can be costly for such individuals to synchronize their activities with others, as they are forced to switch between eating, lying down, or moving at a communal time rather than at a time that is optimal for them as individuals.^{13} Alternatively, a group may split into subgroups that consist of individuals with similar switching times (such as all males and all females, or juveniles and adults), and then the costs of synchronization are lower.^{13–17} Such synchronization costs depend on the different activities of animals in a group, so some animals (e.g., baboons) break up into subgroups for foraging, particularly when food is scarce, and then come together into larger groups for sleeping.^{18}

Social splitting between two categories (e.g., males–females or calves–adults) has been examined using an ordinary-differential-equation (ODE) model, whose performance was tested using data on mixed-sex grouping in red deer.^{12} However, even for animals in the same category (e.g., males), activity synchronization can vary significantly, as it can depend on the age, body mass, and health of animals. Consequently, category-based splitting can lead to groups in which animals are still heterogeneous across many other categories. Splitting of animals in different categories can also be seasonal; for example, in nature, mixed-sex social grouping does not occur during the mating season.^{19}

Communal decisions in herds are made either despotically by a dominant animal (or dominant animals) or democratically by the majority of individuals in a group,^{20–22} and the corresponding groups are called “despotic groups” and “democratic groups,” respectively. Modeling of synchronization costs has suggested that costs for despotic groups tend to be higher than those for democratic groups.^{20}

The rest of our paper is structured as follows. In Sec. II, we discuss biological modeling principles and the construction of a CF, which encapsulates the demands of hunger and lying desire of cattle, for groups of cows to stay together or break apart. In Sec. III, we describe a method for determining the demands of hunger and lying desire using a CF and an evolution scheme (ES) that describes changes in the states (eating, standing, and lying down) of the cows. In Sec. IV, we examine the dynamics of cows and study the CF for various parameter values. In Sec. V, we present two examples: (1) group-switching dynamics of cows when a herd that consists of adults splits into a maximum of three groups, and (2) a scenario in which a herd that consists of an even mixture of males and females splits into a maximum of two groups. In Sec. VI, we discuss our results and present ideas for future work.

## II. BIOLOGICAL MODELING PRINCIPLES

We consider the behavior of cows (*Bos taurus*), which make many daily decisions about staying with or leaving a herd. Cows have a two-stage feeding process that involves first grazing (standing up) and then ruminating (predominantly lying down). Together, lying down and ruminating can occupy up to 65% of a cow's day.^{23,24} Both grazing and lying (including ruminating) are essential for successful digestion of grass,^{25} but cows have to stop these actions if their herd decides to move to another area; this can occur 15–20 times a day.^{24} Each individual cow has similar—but not identical—needs for lying and grazing,^{23,26} so keeping up with a herd each time it moves can be considerably costly because of interrupted grazing or lying times. This cost can include a reduced growth rate in young cattle^{27,28} and physiological and behavioral symptoms of “stress” when a cow is deprived of adequate opportunities for lying down.^{29,30}

Reference 13 examined costs from synchronization, as animals often need to change their behavior (e.g., staying in one place versus moving to another place) at a communal time rather than at their ideal time. In our work, we consider both a synchronization cost and a cost due to predation risk. We assume that large groups encounter a large synchronization cost and small groups increase the cost of predation risk.^{7,8} Therefore, an “optimal” group size is neither too large nor too small. Moreover, we assume that the synchrony can vary within groups, so one set of cows can be eating while another set of cows is lying down or walking (in the neighborhood of others).

We construct a cost function (CF) based on the following four principal assumptions:

Herds are fully democratic when cows take communal decisions, as this reduces cost.

^{20}Cows are free to switch between groups, which freely form or dissolve.

^{31–33}Fission of groups depends only on cows' hunger, lying desire, and predation risk.

The predation risk of a group is an exponential function of the inverse of the group size.

The decrease of predation risk with group size in assumption (iv) arises from the facts that having more animals in a group contributes to greater vigilance,^{1,2,34} a higher dilution effect,^{35,36} and a higher confusion effect.^{37,38} Consequently, a larger group size tends to result in a lower predation risk. Motivated by empirical studies in Refs. 34 and 39–41, which described an inverse exponential relationship between group size and predation risk, we use such a relationship in assumption (iv).

We model the CF, which we denote by *C* in Sec. III B, as a convex combination of costs from hunger (*h*), desire to lie down (*f*), and predation risk (*r*). We thus write

where $\lambda ,\mu \u2208[0,1]$ are parameters. In Eq. (1), “hunger” refers to the grazing demand of cows in a group, and “lying desire” is their demand to lie down. We compute their hunger and lying desire at each time step using a previously-introduced evolution scheme (ES)^{42} for cows to change their state (where eating, standing, and lying down are the three possible states), and we quantify synchronization cost based on cows' hunger and desire to lie down. We assess the cost from hunger (respectively, lying desire) as the mean over all groups of the standard deviation of hunger (respectively, lying desire) within each group, and we model the cost from risk as a function of the group size. During each time step, we minimize the CF over all groups that one can construct from a given herd, where we specify a maximum number of groups, and determine the lowest splitting cost. We determine the optimum group sizes using the minimum of the CF, as it rewards groups with homogeneous demands for hunger and lying desire. This construction enforces perfect synchronization of activity within each group. Our modeling framework is very flexible, and we can examine more general situations by employing different CFs, measuring synchrony in different ways, and considering other extensions.

## III. TEMPORAL EVOLUTION AND MODELING GROUP SPLITTING

As in Ref. 42, when considering a herd, we simulate cows' hunger (i.e., desire to eat) and lying desire (i.e., desire to lie down) and changes of state between eating, lying down, and standing. We then present a CF and optimize it to determine the lowest-cost splitting of the herd.

### A. Temporal evolution and change of states of cows

Cows interact with each other through the ES, which helps provide some understanding of their cooperative activities. We augment the ES in Ref. 42 by formulating it as an iterative scheme that we combine with our CF. In this model, each individual cow is a piecewise-smooth dynamical system, and a cow switches between three discrete states: eating ($E$), lying down ($R$), and standing ($S$). There are also continuous variables, $x\u2208[0,1]$ and $y\u2208[0,1]$, that, respectively, represent the cows' desire to eat and desire to lie down. The dynamics of a single cow are given by state-switching conditions and the following set of differential equations:

where $\xi \u2032i$$is\u2009the\u2009rate\u2009of\u2009increase\u2009of\u2009hunger\u2009,$$\xi i\u2033$$is\u2009the\u2009decay rate of\u2009hunger\u2009,$$\zeta \u2032i$ is the rate of increase of desire to lie down, and $\zeta i\u2033$$is\u2009the\u2009decay\u2009rate\u2009of\u2009desire\u2009to\u2009lie\u2009down$ of the *i*th cow. The parameters $\xi \u2032i,\u2009\xi i\u2033,\u2009\zeta \u2032i$, and $\zeta i\u2033$ are all positive. These parameters can be different for different cows, although for simplicity we did not include the subscript *i* in Eq. (2). If two cows have similar parameter values, we expect them to exhibit similar dynamics. Based on the hypothesis that it is good for cows to eat when other cows are eating and to lie down when other cows are lying down, one can extend the “single-cow model” in Eq. (2) into a coupled dynamical system by allowing the individual cows to interact, and we use a time-dependent adjacency matrix to encode which cows are interacting with each other (see Sec. III C). In Eq. (3) below, we indicate how coupling influences the dynamics of cows.

As we mentioned previously, we modify the coupled system in Ref. 42 to produce an iterative scheme. To simplify our exposition (though at the cost of some technical correctness in the context of animal behavior), we sometimes use the terms “lying desire” to represent “desire to lie down” and “hunger” to represent “desire to eat.” We study the dynamics of the cows at each instance when the state variable changes from one state to another, and we record *x* and *y* for the cows only at these times. Thus, for $t\u2208{1,\u2026,T\u22121}$ and $i\u2208{1,\u2026,n}$, the discrete-time variables $xi(t)\u2208[0,1]$ and $yi(t)\u2208[0,1]$, respectively, represent the level of hunger and desire to lie down of the *i*th cow when the discrete-time state variable $\theta i(t)$ changes at time *t*. The variable $\theta i(t)$ represents the new state of cow *i* at time *t*; it can be eating ($E$), lying down ($R$), or standing ($S$).

As one can see from the paragraph above, the *i*th cow is described by three variables: $\theta i(t),\u2009xi(t)$, and $yi(t)$. For times $t\u2208{1,\u2026,T\u22121}$ and cows $i\u2208{1,\u2026,n}$, the time-dependent coupling of cows is given by the differential equations

where

with

The time-dependent adjacency matrix $A(t)=[aij(t)]n\xd7n$ represents a network of cows at time *t*. Its components are

Thus, $di(t)=\u2211j=1naij(t)$ is the degree (i.e., the number of other cows with which it interacts) of cow *i*. We will discuss such interactions in terms of cow groupings in Sec. III C. The nonnegative parameters *σ _{x}* and

*σ*, respectively, represent coupling strengths with respect to hunger and desire to lie down.

_{y}The switching condition of the state variable $\theta i(t)$ of the *i*th cow at time step *t* is

where we use the parameter $\delta \u2208(0,1)$ to exclude the point $(xi(t),yi(t))=(0,0)$ from the variable domain (because it creates degenerate solutions). Note in Eq. (7) that cow *i* switches to the standing state if *either* its hunger reduces to δ or its desire to lie down reduces to δ (as long as neither of them has a value of 1).

We study the dynamics of cows at discrete times, so we examine a Poincaré section that we construct (using ideas from Ref. 43) by considering switches between different states. We can thereby study the dynamics given by Eqs. (3) and (7). See the schematic in Fig. 1. The boundaries of this Poincaré section are:

These four boundaries arise from the switching conditions in Eq. (7); the first pair of conditions yields the first two boundaries, and the second pair yields the last two boundaries. At time *t*, the variables $xi(t),\u2009yi(t)$, and $\theta i(t)$ of the *i*th cow are represented by one of the boundaries, and then the cow switches to another boundary in the subsequent time step according to the switching condition in Eq. (7).

We solve the dynamical system in Eq. (3) for *n* cows in *T* time steps together with the switching conditions in Eq. (7). The solution gives the discrete dynamics of the *i*th cow in terms of $xi(t),\u2009yi(t)$, and $\theta i(t)$ at each time step *t*. We show the derivation of these solutions in the Appendix as an iterative scheme. As one can see in the left panel of Fig. 1, at time step *t*, each cow is in one of three states ($E,\u2009R$, or $S$) at the start of the time step, and it switches to one of the other two states, where it starts the $(t+1)$th time step. The last two equations in Eq. (8) collectively explain the standing state, so both the lower and the left boundaries of the Poincaré section (see the right panel of Fig. 1 and also Fig. 2) represent the standing state. Thus, in the Poincaré section, the starting point of each cow at a given time step is on one of four boundaries, and the end point at that time step is on a boundary that represents a new state (for which there are two possibilities). We present the corresponding iterative scheme of the solution in Table I, in which we use the following notation:

where $i\u2208{1,\u2026,n}$ and $t\u2208{1,\u2026,T}$.

Case 1: [Eqs. (A1) and (A4)] | If $xi(t)=1,\u2009\delta \u2264yi(t)\u22641$, and $\theta i(t)=E$ |

Subcase a: if $yi(t)\u2265\delta \gamma i\u2032\eta i\u2033$, then $xi(t+1)=[yi(t)]\eta i\u2033\gamma i\u2032,\u2009yi(t+1)=1$, and $\theta i(t+1)=R$ | |

Subcase b: if $yi(t)<\delta \gamma i\u2032\eta i\u2033$, then $xi(t+1)=\delta ,\u2009yi(t+1)=\delta \u2212\gamma i\u2032\eta i\u2033yi(t)$, and $\theta i(t+1)=S$ | |

Case 2: [Eqs. (A2) and (A5)] | If $\delta \u2264xi(t)<1,\u2009yi(t)=1$, and $\theta i(t)=R$ |

Subcase a: if $xi(t)\u2265\delta \eta i\u2032\gamma i\u2033$, then $xi(t+1)=1,\u2009yi(t+1)=[xi(t)]\gamma i\u2033\eta i\u2032$, and $\theta i(t+1)=E$ | |

Subcase b: if $xi(t)<\delta \eta i\u2032\gamma i\u2033$, then $xi(t+1)=\delta \u2212\eta i\u2032\gamma i\u2033xi(t),\u2009yi(t+1)=\delta $, and $\theta i(t+1)=S$ | |

Case 3: [Eqs. (A3) and (A6)] | If $xi(t)=\delta ,\u2009\delta \u2264yi(t)<1$, and $\theta i(t)=S$ |

Subcase a: if $yi(t)\u2264\delta \gamma i\u2032\eta i\u2032$, then $xi(t+1)=1,\u2009yi(t+1)=\delta \u2212\gamma i\u2032\eta i\u2032yi(t)$, and $\theta i(t+1)=E$ | |

Subcase b: if $yi(t)>\delta \gamma i\u2032\eta i\u2032$, then $xi(t+1)=[yi(t)]\u2212\eta i\u2032\gamma i\u2032\delta ,\u2009yi(t+1)=1$, and $\theta i(t+1)=R$ | |

Case 4: [Eqs. (A3) and (A7)] | If $\delta <xi(t)<1,\u2009yi(t)=\delta $, and $\theta i(t)=S$ |

Subcase a: if $xi(t)\u2265\delta \eta i\u2032\gamma i\u2032$, then $xi(t+1)=1,\u2009yi(t+1)=[xi(t)]\u2212\gamma i\u2032\eta i\u2032\delta $, and $\theta i(t+1)=E$ | |

Subcase b: if $xi(t)<\delta \eta i\u2032\gamma i\u2032$, then $xi(t+1)=\delta \u2212\eta i\u2032\gamma i\u2032xi(t),\u2009yi(t+1)=1$, and $\theta i(t+1)=R$ |

Case 1: [Eqs. (A1) and (A4)] | If $xi(t)=1,\u2009\delta \u2264yi(t)\u22641$, and $\theta i(t)=E$ |

Subcase a: if $yi(t)\u2265\delta \gamma i\u2032\eta i\u2033$, then $xi(t+1)=[yi(t)]\eta i\u2033\gamma i\u2032,\u2009yi(t+1)=1$, and $\theta i(t+1)=R$ | |

Subcase b: if $yi(t)<\delta \gamma i\u2032\eta i\u2033$, then $xi(t+1)=\delta ,\u2009yi(t+1)=\delta \u2212\gamma i\u2032\eta i\u2033yi(t)$, and $\theta i(t+1)=S$ | |

Case 2: [Eqs. (A2) and (A5)] | If $\delta \u2264xi(t)<1,\u2009yi(t)=1$, and $\theta i(t)=R$ |

Subcase a: if $xi(t)\u2265\delta \eta i\u2032\gamma i\u2033$, then $xi(t+1)=1,\u2009yi(t+1)=[xi(t)]\gamma i\u2033\eta i\u2032$, and $\theta i(t+1)=E$ | |

Subcase b: if $xi(t)<\delta \eta i\u2032\gamma i\u2033$, then $xi(t+1)=\delta \u2212\eta i\u2032\gamma i\u2033xi(t),\u2009yi(t+1)=\delta $, and $\theta i(t+1)=S$ | |

Case 3: [Eqs. (A3) and (A6)] | If $xi(t)=\delta ,\u2009\delta \u2264yi(t)<1$, and $\theta i(t)=S$ |

Subcase a: if $yi(t)\u2264\delta \gamma i\u2032\eta i\u2032$, then $xi(t+1)=1,\u2009yi(t+1)=\delta \u2212\gamma i\u2032\eta i\u2032yi(t)$, and $\theta i(t+1)=E$ | |

Subcase b: if $yi(t)>\delta \gamma i\u2032\eta i\u2032$, then $xi(t+1)=[yi(t)]\u2212\eta i\u2032\gamma i\u2032\delta ,\u2009yi(t+1)=1$, and $\theta i(t+1)=R$ | |

Case 4: [Eqs. (A3) and (A7)] | If $\delta <xi(t)<1,\u2009yi(t)=\delta $, and $\theta i(t)=S$ |

Subcase a: if $xi(t)\u2265\delta \eta i\u2032\gamma i\u2032$, then $xi(t+1)=1,\u2009yi(t+1)=[xi(t)]\u2212\gamma i\u2032\eta i\u2032\delta $, and $\theta i(t+1)=E$ | |

Subcase b: if $xi(t)<\delta \eta i\u2032\gamma i\u2032$, then $xi(t+1)=\delta \u2212\eta i\u2032\gamma i\u2032xi(t),\u2009yi(t+1)=1$, and $\theta i(t+1)=R$ |

### B. Cost function (CF) determining the splitting of herds

In this section, we determine the lowest-cost grouping of cows by minimizing a CF. This gives the total number of groups and the number of cows in each group. We suppose that a herd of cows splits into a maximum of *L* distinct groups $N1(t),\u2026,NL(t)$, with $|Nl(t)|=nl(t)>0$ cows in the *l*th group (where $l\u2208{1,\u2026,L}$), at each time step $t\u2208T$. If the herd splits into $L1<L$ groups at some time step *t*_{1}, we set $|Nl(t1)|=0$ for $l\u2208{L1+1,\u2026,L}$.

Our CF is the sum of two components: a synchronization component (SC) and a risk component (RC). The synchronization component (SC) models the cost due to variation in the lying desire and/or hunger of cows, and the risk component (RC) models the cost from predation risk.

#### 1. Synchronization component

Recall from Sec. III A that cow *i*'s hunger is $xi(t)\u2208[0,1]$ and lying desire is $yi(t)\u2208[0,1]$. We re-index the variables $xi(t)$ as $xk,l(t)$ and $yi(t)$ as $yk,l(t)$, respectively, to denote the hunger and lying desire of the *k*th cow in the *l*th group at the *t*th time step. Because hunger and lying desire are separate motivations in cows, we compute the two groupings independently, so that cows are optimally homogeneous with respect to hunger (case I) or optimally homogeneous with respect to lying desire (case II). Of these two groupings, we then select the one with the lower synchronization cost.

**Case I**: We sort cows according to increasing hunger, and we place the first $n1(t)$ cows into group $N1(t)$, the next $n2(t)$ cows into group $N2(t)$, and so on. In each group, the synchronization cost from hunger represents the heterogeneity of hunger within the group. As a simple way to quantify this cost, we use the mean of the standard deviations of cows' hunger within the groups and thus calculate

to assess the SC due to hunger. Similarly, we quantify the heterogeneity of groups with respect to lying desire as the mean of the standard deviations of cows' lying desire in groups by calculating

**Case II**: Similar to case I, we sort cows according to increasing lying desire and place them into groups $N1(t),\u2026,NL(t)$. We again compute hunger and lying desire from the means of the standard deviations within the groups, and we denote them by $h2(t)$ and $f2(t)$, respectively.

From the cow groups that we find in cases I and II, we choose the grouping that yields the minimum SC. The hunger $h(t)$ and the lying desire $f(t)$ of the cow herd at time *t* are thus

#### 2. Risk component

Unlike hunger and lying desire, a herd's predation risk is independent of the individuals' states and depends only on the group size. The group size is related inversely to the risk of being attacked by predators,^{1,2,34–38} and we model the predation risk $rl\u2208(0,1]$ of the *l*th group (which has size $nl>0$) as an inverse exponential function of the group size:^{34,39–41}

where *c* is a constant. We assume that the predation risk is small when a group has sufficiently many cows, and we use an associated condition to compute the constant *c*. We denote this sufficient group size (the so-called “safe size”) by *n _{s}*, and we denote the small risk to which the risk function converges (the so-called “safety level”) at this size by

*τ*. The constant

*c*is thus $\u2212(1\u2212ns)/ln(\tau )$, so the RC of the group is

In Fig. 3, we show the relationship given by Eq. (14) for a group with safety level *τ* and safe size *n _{s}*.

We compute the predation risk of each group, and we treat its mean

as the risk of the herd.

In real situations, the safe size and safety level depend on the environment in which a herd lives. If the environment is either dense with predators or vulnerable to predation, the safe size should be comparatively large to achieve a significant safety level. As an example, we use Eq. (15) and compute the risk of splitting a herd of *n* = 20 cows into two groups with a safety level of $\tau =1/30$ and safe sizes of *n _{s}* = 10,

*n*= 20, and

_{s}*n*= 30 (see Fig. 4). We thereby illustrate that large safety sizes model riskier situations for a herd than small safety sizes, independently of how the herd splits. For all safety sizes, we achieve the lowest cost when the herd remains intact (i.e., no splitting), because larger group sizes entail safer herds. We achieve the second-lowest cost when the herd splits into equal-sized groups.

_{s}#### 3. Cost function

We formulate the CF as a convex combination of the costs from hunger, lying desire, and risk of predation:

where $\lambda ,\mu \u2208[0,1]$. For a given herd, which we denote by the set *N*, and a maximum number *L* of groups into which it can split, we minimize (16) over all plausible groups that can be created, and we thereby determine the lowest-cost splitting.

### C. Cost function and temporal evolution

We examine the CF simultaneously with the ES for times $t\u2208{1,\u2026,T}$. At each time step, we update the adjacency matrix $A(t)=[aij(t)]n\xd7n$ in the scheme so that it agrees with the best grouping provided by the optimization of the CF in the previous time step. That is,

The adjacency matrix, which encodes the network architecture of a herd, is an input in Ref. 42. However, in this paper, we update the adjacency matrix at each time step based on an optimum grouping. At each time step, optimizing the CF outputs a lowest-cost grouping until we reach a stopping criterion, which we take to be the maximum time *T*. In Fig. 5, we show a flow chart of this process.

## IV. EXPLORATION OF PARAMETER SPACE

We now explore the effects of parameters on our model. We first examine the dynamics of cows for different coupling strengths, and we then study the CF for different values of the parameters *σ _{x}*,

*σ*,

_{y}*τ*, and

*n*and different coupling strengths.

_{s}### A. Cow dynamics

We explore the dynamics of cows with respect to coupling strength by examining hunger and lying desire on the boundary of a Poincaré section. We then compute the mean group size of a herd for different safety levels.

We undertake these explorations using one herd of *n* = 12 cows that splits into a maximum of three groups (*L* = 3). We simulate hunger and lying desire using the ES (see Sec. III A) followed by computing the CF (16) and optimizing it to determine a lowest-cost grouping at each time $t\u2208T$. As we will discuss shortly, we draw some of the initial conditions and parameter values from probability distributions.

In the ES, we set the initial states of cows to be $\theta i(0)\u2208U{E,R,S}$ for $i\u2208{1,\u2026,n}$, where $U$ denotes a uniform probability distribution over the set in its argument. We add noise sampled from a uniform distribution into the initial conditions and parameters, as it is the simplest type of noise to consider. We determine the initial conditions $xi(0)$ and $yi(0)$ as follows:

For $\theta i(0)=S$, each of the two subcases in Eq. (18) has a 50% chance of being the initial condition. We also make the following parameter choices for the ES: $\xi \u2032i\u2208U[0.0995,0.1005]$, $\xi i\u2033\u2208U[0.1495,0.1505],\u2009\zeta \u2032i\u2208U[0.0495,0.0505],\zeta i\u2033\u2208U[0.1995,0.2005]$, and $\delta =.25$.

Cows are social animals, and their behavior is influenced by what other cows are doing.^{44} We model the strength of such interactions mathematically using the coupling parameters *σ _{x}* and

*σ*in the ES. Using different values for these parameters in different simulations allows us to examine different biological scenarios, such as strongly interacting cows versus weakly interacting cows, and it can be helpful for understanding the dynamics of state changes of cows in these different scenarios. As an initial example, we let $\sigma x=0$ and $\sigma y=0$ (i.e., uncoupled cows) and run the ES for time

_{y}*T*= 400 to simulate hunger $xi(t)$ and lying desire $yi(t)$ for $i\u2208{1,\u2026,12}$ and $t\u2208{1,\u2026,400}$. We also set

*n*= 4 and $\tau =0.2$ in the RC;

_{s}*L*= 3 in the SC; and $\lambda =0.33$ and $\mu =0.33$ in the CF. We then consider three other values of the coupling strengths: $(\sigma x,\sigma y)=(0.05,0.05),\u2009(\sigma x,\sigma y)=(0.2,0.2)$, and $(\sigma x,\sigma y)=(0.6,0.6)$. In each case, we simulate $xi(t)$ and $yi(t)$ for $i\u2208{1,\u2026,12}$ and $t\u2208{1,\u2026,400}$. We determine $xi(0)$ and $yi(0)$ from Eq. (18) with initial states $\theta i(0)\u2208U{E,R,S}$ for $i\u2208{1,\u2026,12}$ a single time (as opposed to determining different values from the same distribution) and perform all four simulations with the same parameter choices (aside from coupling strengths). For each of the four cases above, we simulate one realization of the dynamics. In Fig. 6(a), we show the hunger and lying desire for cow

*i*= 1. The uncoupled case is in the top-left panel, and the coupled cases are in the top-right panel ($(\sigma x,\sigma y)=(0.05,0.05)$), bottom-left panel ($(\sigma x,\sigma y)=(0.2,0.2)$), and bottom-right panel ($(\sigma x,\sigma y)=(0.6,0.6)$).

Unsurprisingly (and by construction), we observe similar dynamics for each cow when they do not interact with each other, but this is not the case when cows are allowed to interact. In other words, $(\sigma x,\sigma y)=(0,0)$ corresponds to modeling cows as independent oscillators, whereas the other cases correspond to coupled oscillators. To compare the dynamics in the four cases $(\sigma x,\sigma y)=(0,0),(\sigma x,\sigma y)=(0.05,0.05),\u2009(\sigma x,\sigma y)=(0.2,0.2)$, and $(\sigma x,\sigma y)=(0.6,0.6)$, we measure the percentage of the length of the Poincaré-section boundary that the first cow's orbit (*x*_{1}, *y*_{1}) intersects in each case. To do this, we discretize each side of the boundary of the Poincaré section into 75 intervals $[0.25,0.26),\u2009[0.26,0.27),\u2009\u2026,\u2009[0.99,1]$, and we then compute the percentage of the number of intervals that the orbit intersects. For $(\sigma x,\sigma y)=(0,0),\u2009(\sigma x,\sigma y)=(0.05,0.05),\u2009(\sigma x,\sigma y)=(0.2,0.2)$, and $(\sigma x,\sigma y)=(0.6,0.6)$, these percentages are about 20.59%, 41.18%, 61.03%, and 66.91%, respectively.

To examine the effect of different coupling strengths in a biological context, we compute the mean group size of the herd (which has *n* = 12 cows) with respect to safety levels and coupling strengths. We first consider $\sigma x=\sigma y=0$ and set the initial states of the ES using $\theta i(0)\u2208U{E,R,S}$ for $i\u2208{1,\u2026,12}$ and the initial conditions using Eq. (18). We use the parameter values $\xi \u2032i\u2208U[0.0995,0.1005]$, $\xi i\u2033\u2208U[0.1495,0.1505],\u2009\zeta \u2032i\u2208U[0.0495,0.0505],\u2009\zeta i\u2033\u2208U[0.1995,0.2005],\u2009\delta =0.25$, and *T* = 400 in the ES; *L* = 3 in the SC; *n _{s}* = 4 in the RC; and $\lambda =0.33$ and $\mu =0.33$ in the CF. We consider 41 safety levels $\tau \u2208{k/40|k=0,\u2026,40}$. For each safety level

*τ*, our simulation generates cow groups ${Nl(t)|l=1,2,3\u2009\u2009and\u2009\u2009t=1,\u2026,400}$. From $|Nj(t)|=nl(t),$ we compute the mean group size and the standard deviation of the group sizes. We perform similar simulations for $\sigma x=\sigma y=0.05,\u2009\sigma x=\sigma y=0.2$, and $\sigma x=\sigma y=0.6$ with the same initial states, initial conditions, and parameter values (other than the coupling strengths). We then compute the mean group size and standard deviation with respect to the safety size in each case.

In Fig. 6(b), we plot the mean group size and group-size standard deviation for each choice of coupling strengths. For each choice, we observe that the mean group size is about 4 for *τ* = 0 and $\tau \u2208[0.8,1]$, that it increases for $\tau \u2208(0,0.58]$, and that it decreases for $\tau \u2208(0.6,0.8)$. We observe [see Eq. (15)] that the herd can maintain a small risk even when splitting into small groups for small safety levels. For larger safety levels, the mean group size must be larger to ensure that the cost of the RC is sufficiently small. However, beyond some value of the mean group size, the SC starts to dominate the CF. The mean group size thus starts to become smaller for larger safety levels. In this example, the value of the safety level at which this trade-off balances is about $\tau =0.6$. We also observe that stronger coupling entails larger mean group sizes. For large safety levels, cows can form large groups with similar dynamics without the herd incurring a significant cost.

### B. Cost function

We now examine how the CF changes with respect to the coupling strengths *σ _{x}* and

*σ*, the safe size

_{y}*n*, and the safety level

_{s}*τ*. We perform three numerical experiments: one to examine the effect of coupling strengths; another to examine the effects of the safe size and safety levels, and another to compare the effect of safety level for zero and nonzero coupling strengths.

We have already seen in Sec. IV A that different coupling strengths yield different dynamics for state switches in cows. From a biological perspective, a larger coupling strength corresponds to stronger interactions between cows, and we wish to explore how different interaction strengths affect the cost of synchronizing behavior. In our simulations, we average the cost over five realizations of parameter values of herds of *n* = 15 cows. Specifically, we generate five sets of initial states using $\theta i(0)\u2208U{E,R,S}$ for $i\u2208{1,\u2026,15}$ and then use Eq. (18) to generate five initial conditions. We then generate five sets of values for the parameters $\xi \u2032i\u2208U[0.0995,0.1005],\u2009\xi i\u2033\u2208U[0.1495,0.1505],\zeta \u2032i\u2208U[0.0495,0.0505]$, and $\zeta i\u2033\u2208U[0.1995,0.2005]$. We set $\delta =0.25,\u2009n=15$, and *T* = 20 in the ES and consider the coupling strengths $\sigma x=\sigma y=k/150$, where $k\u2208{1,\u2026,30}$. We also set *L* = 3 in the SC; *n _{s}* = 4 and $\tau =0.2$ in the RC; and $\lambda =0.33$ and $\mu =0.33$ in the CF. We then run the dynamics for each initial condition and compute five cost values for each choice of coupling strengths $\sigma x=\sigma y=k/150$. In Fig. 7(a), we plot the standard deviations of the costs, and we observe that they decrease with

*σ*(and hence with

_{x}*σ*) until appearing to saturate once the coupling strength reaches a value of about 0.065.

_{y}Our model assesses the effect of risk in the cost using the RC. The risk levels in the RC depend on both the safety level *τ* and the safe size *n _{s}* [see Eq. (15)]. In risky environments, we expect the values

*τ*and

*n*to be larger than in safe environments. To assess the influence of these parameters on the CF, we perform simulations with the initial states $\theta i(0)\u2208U{E,R,S}$ for $i\u2208{1,\u2026,15}$ and determine the other initial conditions from Eq. (18). We use $\xi \u2032i\u2208U[0.0995,0.1005],\u2009\xi i\u2033\u2208U[0.1495,0.1505],\u2009\zeta \u2032i\u2208U[0.0495,0.0505],\u2009\zeta i\u2033\u2208U[0.1995,0.2005],\u2009\delta =0.25$, $\sigma x=0.1,\u2009\sigma y=0.1$,

_{s}*n*= 20, and

*T*= 20. Additionally, we let

*L*= 3 in the SC; $ns=k1$ (with $k1\u2208{2,\u2026,20}$) and $\tau =0.05k2$ (with $k2\u2208{1,\u2026,19}$) in the RC; and $\lambda =0.33$ and $\mu =0.33$ in the CF. For each $(k1,k2)\u2208{2,\u2026,20}\xd7{1,\u2026,19}$, we perform one realization (so we only consider one value of each parameter determined from a probability distribution). For each

*τ*and

*n*, we run the simulation for

_{s}*T*= 20 times steps and consider the cost at each time. In Fig. 7(b), we show the mean cost as a function of safe size

*n*and safety level

_{s}*τ*. We observe that the cost is low for small values of

*n*and

_{s}*τ*, and that it becomes progressively larger for larger parameter values until it appears to saturate.

In our two experiments above, we examined how the CF depends on coupling strength and safety level. We now examine the temporal variation of the CF versus the safety level for both uncoupled cows and coupled cows. We generate one set of initial states using $\theta i(0)\u2208U{E,R,S}$ for $i\u2208{1,\u2026,15}$ and one set of initial values using Eq. (18). We then choose the parameters in the ES, SC, and CF as in our simulations above to examine the influence of coupling strengths on the CF. We set $\sigma x=0$ and $\sigma y=0$ for the uncoupled cows and $\sigma x=0.1$ and $\sigma y=0.1$ to examine a situation with coupled cows. In the RC, we let *n _{s}* = 4 and consider a safety level of $\tau =k/20$, where $k=0,\u2026,20$. In Fig. 7(c), we show the cost as a function of time and safety level for both uncoupled and coupled cows. We observe for uncoupled cows that the cost is larger for a larger safety level. Importantly, however, this need not be the case for coupled cows.

## V. BIOLOGICALLY-MOTIVATED EXAMPLES

We examine the CF (16) using two biological examples: (1) a herd that splits into up to three groups and (2) a herd with males and females that splits into two groups.

### A. Example 1

In this example, we illustrate a scenario of a herd splitting into up to three groups. It also helps convey the effect of choosing parameter values in Eq. (16) and the relationship between groupings and their associated costs.

We consider a herd of *n* = 12 cows that we allow to split into a maximum of *L* = 3 groups during *T* = 30 time steps. We first simulate hunger and lying desire, then compute the CF, and finally optimize the CF to determine the lowest-cost grouping at each time. We consider a single realization of the model (i.e., one example herd) and use it to illustrate the general notion of trade-offs in the CF.

In the ES, we set the initial states of the cows to be $\theta i(0)\u2208U{E,R,S}$ for $i\u2208{1,\u2026,n}$, and we recall that $U$ denotes a uniform probability distribution over the set in its argument. We set the initial conditions $xi(0)$ and $yi(0)$ according to Eq. (18). We also make the following parameter choices for the ES: $\xi \u2032i,\u2208U[0.0995,0.1005],\u2009\xi i\u2033\u2208U[0.0495,0.0505],\u2009\zeta \u2032i\u2208U[0.1245,0.1255]$, $\zeta i\u2033\u2208U[0.0745,0.0755],\u2009\delta =0.25,\sigma x=0.1$, and $\sigma y=0.1$. We set the parameters in the CF and RC to be *n _{s}* = 4, $\tau =0.2,\u2009\lambda =0.2$, and $\mu =0.2$.

A herd of 12 cows can split into a maximum of 3 groups in 19 different combinations of group sizes (see Table II). We assign an index for each combination to simplify the labeling in our figures. We also run the ES together with the CF for another two instances of the CF parameters: $\lambda =0.6,\u2009\mu =0.2$ and $\lambda =0.2,\u2009\mu =0.6$. We show our results at time *t* = 20 for all three examples in Fig. 8. In the figure, the highest risk occurs for $n1,2,3(20)=10,1,1$, in which two individual cows have separated from a herd. The second-highest risk occurs when $n1,2,3(20)=0,1,11$, in which one cow has separated from a herd. The lowest risk occurs when the entire herd stays together (index 1) or when it splits into equal groups (index 7), where we note that the group size of 6 is larger than the safety size *n _{s}* = 4. One can consider equally-weighted cost components in the convex combination that constitutes the CF or change the importance of components by increasing the weight of hunger [see Fig. 8(a)], lying desire [see Fig. 8(b)], or risk [see Fig. 8(c)].

Index . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | 9 . | 10 . | 11 . | 12 . | 13 . | 14 . | 15 . | 16 . | 17 . | 18 . | 19 . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

$n1(t)$ | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 10 | 9 | 8 | 7 | 6 | 8 | 7 | 6 | 5 | 6 | 5 | 4 |

$n2(t)$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 1 | 2 | 3 | 4 | 5 | 2 | 3 | 4 | 5 | 3 | 4 | 4 |

$n3(t)$ | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | 3 | 3 | 4 |

Index . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | 9 . | 10 . | 11 . | 12 . | 13 . | 14 . | 15 . | 16 . | 17 . | 18 . | 19 . |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

$n1(t)$ | 12 | 11 | 10 | 9 | 8 | 7 | 6 | 10 | 9 | 8 | 7 | 6 | 8 | 7 | 6 | 5 | 6 | 5 | 4 |

$n2(t)$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 1 | 2 | 3 | 4 | 5 | 2 | 3 | 4 | 5 | 3 | 4 | 4 |

$n3(t)$ | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | 3 | 3 | 4 |

We now examine the temporal grouping in this scenario with parameter values $\lambda =0.33$ and $\mu =0.33$. In Fig. 9(a), showing 6 arbitrary cows out of 12 in total, we see that cows freely switch their groups to achieve the optimum value of the CF (16). The cow that we represent with purple crosses switches between two groups during the simulation, whereas the other five cows switch between all three groups. In Fig. 9(b), we show the total number of groups in the herd, which consists of a single group at times *t* = 19 and *t* = 23 and consists of three groups at times *t* = 3, *t* = 7, *t* = 21, *t* = 24, *t* = 25, and *t* = 28. In Fig. 9(c), we show the total cost and thereby reveal that it can be more costly for the herd to stay together as a single group than to split up (at times *t* = 19 and *t* = 23). We also note the low costs for times *t* = 3, *t* = 7, *t* = 21, *t* = 25, and *t* = 28, when the herd consists of three groups. Note that we have illustrated trade-offs in the CF specifically for the initial condition and parameter values in our example, and we expect to see qualitatively different trade-offs for different initial conditions and parameter values. (Additionally, the “high” and “low” costs are not much different from each other.) However, the notion of such trade-offs is a rather general one.

### B. Example 2

We now examine mixed-sex grouping dynamics in a herd that consists of two distinct categories of adult cows: males and females. This type of grouping is known to occur in some animal groups (e.g., in red deer^{12}), so we study the same phenomenon in our model of cow herds. Adult male cows require more energy and rest than female cows, as the former tend to have larger body weights.^{45,46} We therefore assume that the males' rates of change of hunger and lying desire are larger than those of females. Mathematically, we implement this assumption by using larger values of the parameters $\xi \u2032i,\xi i\u2033,\zeta \u2032i$, and $\zeta i\u2033$ of cows in the male group than for those in the female group. (It is reasonable that, e.g., a male cow becomes hungrier faster than a female cow, but our analogous assumption is much less reasonable for the sating of hunger and the desire to lie down.)

We consider a herd of 10 cows of two different types. There are five cows (where $i\u2208{1,\u2026,5}$ indexes the cow) with large body weights, and the remaining five cows ($i\u2208{6,\u2026,10}$) have small body weights. As in Sec. V A, we simulate the hunger and lying desire of cows with the ES (see Sec. III A) and determine a lowest-cost grouping by optimizing the CF (16). We set the initial states of cows of the first and second types as eating and lying down, respectively. For cows of a given type, the variables have very similar initial values. Specifically, they are the same, except that we perturb them additively with a small amount of uniform noise

where $\varphi i,\varphi \u2032i\u220810\u22123U[0,1]$ and $\delta =0.25$. We choose uniform additive noise because it is the simplest type of noise to consider. We set $\sigma x=0.2$ and $\sigma y=0.2$ in the ES, and we determine the other parameters so that the first group consists of cows with a large body mass and the second group consists of cows with a small body mass. That is,

We set the parameters in the CF and RC to be *n _{s}* = 3, $\tau =0.2,\u2009\lambda =0.33$, and $\mu =0.33$. We run the ES for

*T*= 30 time steps, and we consider the value of the CF at each step. As in our examples in Sec. V A, we use only one realization, and we note that the noise in Eq. (19) has a small magnitude. During time steps 0–10 and 20–28, we see in Fig. 10(a) that all of the cows are in groups with the other cows of their own sex (i.e., with others of similar sizes, hunger, and desire to lie down). However, during time steps 11–19 and 29–30, some cows are not in their “proper” group, and the cost becomes high [see Fig. 10(c)], although the CF minimizes the cost to achieve a lowest-cost grouping. We show the number of mismatched cows in the groups in Fig. 10(b). We observe that the cost is large when cows are in mismatched groups, but it is low when cows are in their proper (i.e., single-sex) groups.

## VI. CONCLUSIONS AND DISCUSSION

We developed a framework for modeling the lowest-cost splitting of a herd of cows by optimizing a cost function (CF) that quantifies their hunger, desire to lie down, and predation risk. Lying in groups offers protection from predators,^{19,47–49} but synchronization can also be costly to individuals, as some portion of a herd has to change behavior to eat or lie down at a communal time rather than at an optimally beneficial time.^{11–13} In this paper, we examined situations in which cow herds split into groups such that cows' hunger and lying desire are relatively homogeneous within a group, while ensuring that further splitting does not result in overly small groups, which would be more vulnerable to predation.

We employed the evolution scheme (ES) from Sun *et al.*^{42} and input cows' time-dependent interactions in terms of a adjacency matrix $A(t)$ that encodes the lowest-cost grouping obtained by optimizing the CF. The adjacency matrix provides an interface between the CF and ES, and our framework can be used with arbitrarily intricate CFs, ESs, and interaction patterns. In Ref. 42, the network architecture $A(t)$, which indicates which cows interact with each other at each time *t*, was imposed as part of the model. In the present paper, however, we took a different approach: we determined $A(t)$ based on an optimal grouping at the previous time step (after imposing a group structure at *t* = 0 as a part of the initial conditions). Because hunger and lying desire are two separate motivations of a cow, we optimized the CF independently for each of them in each time step to obtain two different groupings, and we then used the grouping with the lower total cost among the two possibilities. For convenience, we imposed a maximum number of groups into which a herd can split, as it reduces the computational complexity of our approach. We assessed the cost contributions from hunger and lying desire using the standard deviation of the associated individual preferences in each group (although one can replace the standard deviation by any measure of dispersion).

In Sec. IV, we first examined how cow dynamics are affected by coupling strengths, and we then examined the CF for different parameter values. We simulated hunger and lying desire of cows for four sets of coupling strengths and observed different dynamics [see Fig. 6(a)] in the four situations. Setting the coupling strengths to 0 implies that each cow behaves independently [see Eq. (3)], so each cow acts as an independent oscillator. In contrast, for positive coupling strengths, cows interact with each other, and a cow herd is then a set of coupled oscillators. To examine the different dynamics from different coupling strengths in a biologically-motivated context, we computed the mean group size versus the safety level for different coupling strengths. We observed in Fig. 6(b) that large coupling strengths permit large groups that consist of cows with similar needs. We also observed that the mean group size of cows first becomes larger for progressively larger safety levels but then becomes smaller after some value of the safety level [see Fig. 6(b)]. Recall that group sizes in a herd also increase with the safety level [see Eq. (3)]. In sufficiently large groups, the synchronization cost starts to dominate the CF for sufficiently large safety levels, and minimizing the CF starts to encourage smaller groups to minimize the cost. Thereafter, the mean group size decreases with the safety level.

We then studied the influence of coupling strengths *σ _{x}* and

*σ*, safe size

_{y}*n*, and safety level

_{s}*τ*on the CF. We observed [see Fig. 7(a)] that the total cost becomes smaller for progressively larger coupling strengths before saturating. In Fig. 7(b), we illustrated that setting the safe size and safety level to low values entails a low cost. Such low parameter values allow cows to gather into small groups of similar cows without incurring a significant risk to a herd. When the cows are uncoupled, the cost increases monotonically with the increasing safety level, but the cost varies non-monotonically with increasing safety levels for coupled cows [see Fig. 7(c)].

In a biologically-motivated example, we examined group fission and the dynamics of cows switching between groups. In that example, we set the initial states of cows arbitrarily, but one can also choose initial states to examine specific scenarios. To consider a relatively homogeneous herd, we used similar parameter values for different individuals, and we observed the dynamics that result from small differences in these parameter values. We considered a single realization of the model, and other initial conditions and parameter values yield different specific trade-offs while illustrating the same essential idea. Our primary hypothesis, that synchronization can be costly, is illustrated by Figs. 9(b) and 9(c). Specifically, synchronization is very costly when the groups are large and heterogeneous. One can explore trade-offs further by considering risk and synchronization costs with different rates of increase with group size.

One can customize the ES by changing the parameters for the rates of increase in hunger or desire to lie down. This versatility allowed us to model a scenario of mixed-sex grouping in a herd. Adult male cows generally possess larger body masses and require more energy and lying time than adult female cows. We implemented this asymmetry among individuals by imposing larger values of the salient parameters for males than for females.^{50} At times, the heterogeneity in motivations for eating and lying down caused the optimal groups to consist of cow groups other than the single-sex groups [see Fig. 10(c)], but usually optimization of the CF yielded single-sex groups. Single-sex grouping occurs commonly in ungulates^{51} (e.g., cows, deer, and sheep) and are especially pronounced in species with large body-size differences between males and females.^{11,12,14,15} In our exploration of sex grouping, we added uniform noise to the initial conditions and parameters, as it is the simplest type of noise to consider.

One can adjust the CF so that it can be used for herding situations in different environments. A safe environment allows small groups in a herd, in contrast to an unsafe environment, which requires large groups to defend themselves against attacks. Our CF imitates a safe environment if the safe size *n _{s}* is large and the safety level is small. One can control the influence of the cost components (hunger, lying desire, and predation risk) on the CF by tuning parameters, and our approach thereby makes it possible to explore different grouping scenarios, such as analyzing the influence of one or more cost components versus the others for group splitting. Our approach is also very flexible, and one can generalize our CF, the ES, and the interactions among animals (through a time-dependent adjacency matrix) to examine a wide variety of scenarios.

In our paper, we determined group size and splitting by optimizing a CF at each time step. However, because optimally-sized groups are not necessarily stable, it is important to explore the idea of introducing a learning process in which one keeps track of optimal group sizes during past time steps. In the present paper, we imposed a maximum number *L* of groups into which a herd can split. In our examples, the value of *L* was either obvious, as in the sex-grouping example, where we used *L* = 2 (males and females), or hypothetical, as in our example with *L* = 3. However, instead of imposing a maximum number of groups in advance, it is also desirable to examine situations in which the number of groups is an unconstrained output to better reveal an optimal number of groups in herd splitting.

In summary, we developed a versatile model of lowest-cost splitting of a herd of animals that allows numerous generalizations in a straightforward way. We illustrated our model by exploring several plausible scenarios, and we believe that our approach has the potential to shed considerable insight on grouping behavior in animals in a wide variety of situations.

## ACKNOWLEDGMENTS

We thank Jie Sun for his valuable comments on this work. E.M.B. and K.G. were supported by the National Science Foundation (DMS-0404778), and E.M.B. was also supported by the Office of Naval Research (N00014–15-1–2093) and the Army Research Office (N68164-EG and W911NF-12–1- 0276).

### APPENDIX: DERIVATION OF THE DISCRETE DYNAMICS ON THE POINCARÉ SECTION

We solve the differential equations in Eq. (8) using the boundary conditions in Eq. (7). For convenience, we substitute Eq. (9) into these differential equations and expand as follows:

when $\theta i(t)=E$_{,}

when $\theta i(t)=R$_{,}

when $\theta i(t)=S$,

We then solve the differential equations in Eqs. (A1)–(A3) on the boundaries $\u2202E,\u2009\u2202R,\u2009\u2202Sx$, and $\u2202Sy$ given by Eq. (8) as follows (where the subscripts of time *t* indicate what state change is occurring):

when $\theta i(t)=E$ (i.e., on $\u2202E$ of the Poincaré section),

when $\theta i(t)=R$ (i.e., on $\u2202R$ of the Poincaré section),

when $\theta i(t)=Sy$ (i.e., on $\u2202Sy$ of the Poincaré section),

when $\theta i(t)=Sx$ (i.e., on $\u2202Sx$ of the Poincaré section),

## References

As we discussed previously, it is reasonable to assume that males experience faster increases in hunger and desire to lie down than females. However, our parameter choices also entail the seemingly unrealistic situation that they also become sated faster.