Evolution of Specialization in Dynamic Fluids

Previously we found mechanical factors involving diffusion and fluid shear promote evolution of social behavior in microbial populations Uppal and Vural (2018). Here, we extend this model to study the evolution of specialization using realistic physical simulations of bacteria that secrete two public goods in a dynamic fluid. Through this first principles approach, we find physical factors such as diffusion, flow patterns, and decay rates are as influential as fitness economics in governing the evolution of community structure, to the extent that when mechanical factors are taken into account, (1) Generalist communities can resist becoming specialists, despite the invasion fitness of specialization (2) Generalist and specialists can both resist cheaters despite the invasion fitness of free-riding. (3) Multiple community structures can coexist despite the opposing force of competitive exclusion. Our results emphasize the role of spatial assortment and physical forces on niche partitioning and the evolution of diverse community structures.

; May and Seger (1986), all fulfill specialized roles that 24 are vital for the functioning of the bigger whole. The evolution of specialization also gives rise 25 to cooperative metabolic interdependencies in microbial populations and can serve as a strong 26 mechanism for community assembly Zelezniak et al. (2015). 27 There are two classes of evolutionary forces moving a population from having one type of 28 individual performing multiple functions -generalism-, towards one that has multiple types of 29 individuals performing distinct functions -specialism-. The first is "incompatible optimas" Solari Specialists (large red and yellow circles) secrete only one of the two public goods. Cheaters (dark green) secrete none of the public goods. All microbes secrete a metabolic waste (small blue circles). (b) Fitness contour plots and types of stable groups in each fitness variant. In the top row, we plot the fitness contours as a function of public good concentrations 1 and 2 . In the bottom row, we show the types of stable groups in each fitness form. The red line represents the contour corresponding to zero fitness. In the AND case, the red line can never cross the 1 , 2 axes, the fitness is negative when either chemical is not present. Therefore the only types of stable groups are generalists or mixed specialists, as shown in the bottom-left panel. In the OR case, the contours form straight lines. A decrease in one chemical is equally compensated by an increase in the other chemical. The zero contour always crosses the 1 , 2 axes. It is therefore also possible to have pure specialist groups in the OR case, as shown in the bottom-right panel. (c) Evolutionary paths between group types. Groups typically move towards less secretion since cheaters have higher invasion fitness than specialists, who, being "half-cheaters", have higher invasion fitness than generalists.
Every microbe secretes a waste molecule that curbs the growth of those nearby. (3) The secretions 94 and bacteria obey the physical laws of fluid dynamics and diffusion. (4) Whether a microbe secretes 95 both, one, or none of the public goods is hereditary, except for random mutations. However, every 96 phenotype emits waste. 97 We study two models separately. (5) In one, which we call AND, access to both kinds of goods is 98 necessary. In the other, which we call OR, both goods contribute to fitness, but the lack of one can 99 be compensated with the other. 100

101
Our work consists of (1) discrete, stochastic agent based simulations and (2) related continuous 102 deterministic equations. In addition, to gain better analytical understanding, we construct (3) a 103 simple effective model that captures the essential details of (1) and (2). 104 105 We construct equations governing the number density of three phenotypes 1 ( , ), 2 ( , ), 3 ( , ) two chemical secretions that are public goods 1 ( , ), 2 ( , ), and a waste compound 3 ( , ), as a function of space and time . 3 ( , ) is the number density of microbes that secrete both kinds of public goods, to which we refer as "generalists". The microbes that secrete only public good number one or two, are denoted by 1 ( , ) and 2 ( , ), to which we refer as "specialists". Those that secrete no public goods are denoted by 0 ( , ), to which we refer as "cheaters".

Continuous Deterministic Equations
Here indices , = 0, 1, 2, 3 label phenotypes, whereas the index = 1, 2, 3 labels chemicals, i.e. the 106 two public goods and waste. Thus, Equation 1 and Equation 2 comprise 7 coupled spatiotemporal 107 equations. 108 In both equations, the first two terms describe diffusion and advection. The flow field ( , ) is a vector valued function of space and time, and includes all information pertaining the flow patterns in the environment. In general, it is obtained by solving separate fluid dynamics equations. Mutations between phenotypes and secretions of each phenotype is governed by two matrices The secretion rate of chemical by phenotype is given by the matrix element , and its decay 109 rate by . The mutation rate from phenotype to is given by . The diagonal elements 110 indicate the rate at which phenotype mutates into a different category. 111 In our model, the secretion of public goods is binary, i.e. a good is either secreted or not. 112 Mutations toggle on and off with probability whether an individual secretes either public good. A 113 mutation can cause a generalist to become a specialist, but two mutations, one for each secretion 114 function, are required for a generalist to become a cheater. Same with back mutations. 115 The fitness function ( ) determines the growth rate of phenotype . We consider two cases separately: when both public goods are necessary for survival (AND) and when the public goods can substitute for one other (OR).
(AND) = 12 1 2 As we see, in both cases, growth rate increases with the local concentration of public goods, 1 , 2 116 and decreases with the concentration of waste, 3 . is the cost of secreting public good , so that 117 growth of phenotype is curbed by an amount proportional to its public good secretion. Note that 118 waste is produced without any cost.

155
Cooperative groups as Turing patterns 156 Through numerical simulations and analytical formulas, we see that the system gives rise to spatially 157 segregated cooperating groups in a certain parameter range, as shown in Figure 2. Spots or stripes 158 in reaction diffusion systems are known as Turing patterns, which form whenever an inhibiting 159 agent diffuses faster than an activating agent. In our model the inhibiting and activating agents are 160 the waste and the public goods. 161 In general, the structure and size of these cooperating groups will vary with physical parameters. 162 We show in Figure  In our simulations, we observe that cooperative groups of microbes, i.e. spots and stripes, grow 169 and fragment, thereby giving rise to new structures of the same type. These spatial structure of 170 these patterns differ between generalists and specialists, and therefore have a strong effect on the 171 evolutionary trajectory of the system.

172
Effects of secretion cost on specialization 173 We next determine the roles of secretion cost , and flow shear on group structure and hence 174 specialization. To see the effect of trade-offs on specialization, we varied the cost of public good 175 secretion and determined when specialization occurs in both AND and OR fitness forms. To simplify 176 our analysis, we set 1 = 2 = . In order for both types of specialists to then coexist, we also set 177 1 = 2 = . Therefore, generalists pay an overall cost of 2 , specialists pay , and cheaters pay 178 no cost. As such, a specialist mutant will invade a generalist group, and a cheater mutant will 179 invade a specialist group. In the absence of spatial structure and flow, the entire population will be 180 dominated by cheaters and will go extinct. 181 What can we say about the competition between different group types (as opposed to between 182 different strains within a group)? Since the fitness benefit per cost diminishes with larger secretion, 183 one might expect that increasing the cost of the goods would favor the specialists over generalists. 184 Counterintuitively, we find the opposite. Specialist groups indeed grow faster and form larger, 185 expansive, and denser groups, which however are at once taken over by mutants. In contrast, 186 generalists form smaller, sparser, weaker groups that fragment often, which limits the spread of 187 mutants. Therefore, at higher cost , the "weak" generalists are able to coexist and even dominate 188 "strong" specialists (Figure 3a,b). 189 In general, a large uniform population is more susceptible to invading mutants. In contrast, 190 when the population is organized as fragmenting patches, the community structure will prevail 191 as long as the fragmentation rate is larger than the invasive mutation rate. Thus, the type, size, 192 growth and fragmentation of the groups ultimately dictates whether generalism, specialism, or a 193 coexistence of group types are evolutionarily stable.

194
Effect of flow patterns on specialization 195 A shearing fluid flow has been shown to modify social behavior by enhancing the group size and 196 fragmentation rate Uppal and Vural (2018). We therefore expect the shear rate to effect the types 197 and sizes of groups that we see at steady state. 198 For constant shear we used a planar Couette flow, with velocity profile and shear rate given as, where max is the maximum flow rate and is the height of the domain. We used periodic boundary The effect of shear is in general non-trivial and will depend on the group structure observed. 202 We find that a shearing flow enhances group fragmentation rate when a species forms distinct spot 203 patterns and enlarges aggregates when they are closer to forming stripes. 204 In Figure 3c, state to a specialist state (Figure 3-video 2). Thus, fluid shear promotes specialization. 208 Since advective fluid flow is something that one can tune in an experimental or industrial setting, 209 it is exciting to think of possibilities where flow is used to control the social evolution of a microbial 210 community. Furthermore, since shear is in general spatially dependent, we can use different velocity 211 profiles to localize this control to different regions.

212
Combined effects of public good benefit, cooperation cost, and competition on 213 evolution of specialization 214 We study how varying public good benefit, production cost, and waste diffusion affect the stability 215 of different community structures (Figure 4). We find in general, that higher waste diffusion and 216 higher public good benefit helps specialists and higher secretion cost favors generalists. The pie 217 charts also show what conditions leads to coexistence of different group types. 218 If waste diffusion is large, self-competition is lower, and specialists can form denser groups 219 without over-polluting themselves. They can then better utilize public goods secreted by their 220 neighbors. If the public good benefit, 12 , is large, specialists also do better since secreting fewer 221 public goods still gives a large benefit (Figure 4-video 1, Figure 4-video 2). 222 As we have already seen, specialization actually occurs more often when trade-offs are small, i.e. 223 at smaller . At higher , generalists are able to coexist with specialists and constitute the majority 224 of the population (Figure 4-video 3, Figure 3-video 1). 225 We also see that can cheaters persist stably with the population when their invasion fitness 226 is lower than the growth rate of producers. This occurs in regions where producers do not form 227 groups but grow either as stripes or homogeneously in space, which happens when public good 228 benefit is large and when secretion costs are low. In this case, cheaters "chase after" producers, 229 which grow into free space (Figure 4-video 4). High waste diffusion also helps cheaters, since 230 they are able to chase producers without over-polluting themselves or their hosts. When their 231 invasion fitness is about equal to the producer growth rate, cheaters take over fully, driving the 232 Figure 3. Effects of cooperation cost and fluid shear on specialization. Points with error bars represent numerical results averaged over 5 runs at simulation time = 2 × 10 6 s in a domain of size 20 mm × 20 mm.
Error bars correspond to one standard deviation from the mean. Solid lines are from our effective theoretical model given in Appendix 2. (a,b) In each fitness variant, we see that specialization and cheating is more abundant for low costs. At lower costs, generalist groups are "too fit" and form large aggregate structures that are more susceptible to mutations. At higher costs, we see generalists out-compete specialists, in the AND case (a), and coexist with specialists in the OR case (b). At low costs, specialists again dominate. (c,d) Effects of fluid shear. A shearing flow causes groups to fragment quicker and stripes to elongate and grow larger. (c) In the AND case, with secretion cost = 0.12, we observe that shear transitions the system from a coexisting state to a specialist state. Here, shear causes generalists groups to elongate and become more susceptible to mutations and causes specialist groups to fragment quicker than generalists. (d) In the OR case, at cost = 0.17, shear again causes a transition from a coexisting state to a pure specialist state. We therefore see in both cases that flow shear will promote specialization. Parameters not shown in plots are given in Table 1.   Pie charts give the population ratios of generalists, specialists, and cheaters for varying secretion cost , waste diffusion , and public good benefit 12 . Population values were obtained by taking a time average over 1 run for each parameter value. The time average taken over time steps = 1 × 10 6 s to = 2 × 10 6 s. When waste diffusion is roughly larger than the public good diffusion, the population can form spatial structures. Under certain conditions, we see coexistence amongst all three -generalists, specialists, and cheaters. At medium costs, cheaters are able to take over faster than groups reproduce, leading to extinction, seen as empty regions. At higher costs, specialists form smaller groups that fragment quicker than cheaters take over and are stable at steady state. Specialists do better overall when waste diffusion is large, since they can then form denser groups without over-polluting themselves. Higher public good benefit, 12 , also helps specialization, since secreting fewer public goods still gives a large benefit. Interestingly, higher benefit also leads to more extinct states, since cheaters can take over quicker. These trends apply to both AND and OR cases. Parameter values for 12 , , and are shown in the plots. Public good diffusion for the AND case is 12 = 20 × 10 −6 cm 2 s −1 , for the OR case 12 = 25 × 10 −6 cm 2 s −1 ; the rest of the parameters are as given in Table 1.     236 We see two regions of extinction: when public good benefit and waste diffusion are large, at 237 medium costs; and when public good benefit and waste diffusion are low, at high costs. The first 238 case is due to cheaters taking over groups, leading to the tragedy of the commons. Interestingly, 239 this occurs more with higher public good benefit. The population of producers becomes "too fit" 240 and more vulnerable to cheating mutations. For the second case, since costs are high and benefits 241 are low, microbes need to form dense groups to utilize enough goods to be stable. However, due to 242 the low waste diffusion, these groups over-pollute themselves and are no longer stable. 243 We see similar trends for both the AND and OR cases. The main distinction between the two 244 being, for the OR case, we predominately see pure specialist groups and only have mixed specialists 245 in the AND case. We do not see many mixed specialists in the OR case since mutations take over 246 generalists groups quicker and stabilize as pure groups, whereas in the AND case, pure groups 247 would die out unless the complementary specialist also evolves in the same group. The AND 248 structure is therefore essential to have true division of labor, where each type of specialist exists 249 equally in the group. 251 We next study the evolution of specialization in Hagen-Poiseuille and Rankine vortex flows. Again, we set the cost parameter to a value where shear makes the biggest difference. As with the case with constant shear (Figure 3c,d), we set for AND fitness, = 0.12 and for OR fitness, = 0.17. For a Hagen-Poiseuille flow in a two-dimensional pipe, the flow rate and shear rate are given by,

Localization of specialization and coexistence in pipe and vortex flows
From our results with a constant shear (Figure 3c,d), we expect higher shear regions of the pipe to be 252 occupied by specialists and lower shear regions to be occupied by generalists. However, we see the 253 opposite to occur (Figure 5). This is due to boundary and second order effects. Generalist groups on 254 the boundary fragment more often, and longer groups are formed in regions of intermediate shear 255 (Figure 5-video 1). The fragmenting generalists groups act as a source for specialists groups in the 256 intermediate regions of the pipe. Near the center of the pipe where the shear rate is low, groups do 257 not fragment as quickly and are taken over by cheaters. We therefore see a coexistence of species 258 across the pipe, with generalists at the boundary, followed by specialists in the intermediate regions, 259 and no one at the center (Figure 5). 260 Next we study evolution in a Rankine vortex. The flow and shear profiles for a Rankine vortex with radius and circulation Γ are given by, The distribution of specialists and generalists in the vortex agrees better with previous results from 261 constant shear (Figure 5). We see generalists persist in regions of low shear and specialists mainly 262 reside in an annular region where shear is large (supplementary video Figure 5-video 2). 263 In either case we see coexistence of different species across the full domain. The local shear 264 profile dictates which species is stable in each region. A varying shear profile can therefore allow 265 for the coexistence of species across an environment. 266 Figure 5. Coexistence of specialists and generalists in pipes and vortices. The local shear profile dictates which species is stable. When shear is spatially varying, we can get coexistence of generalists and specialists. (a) For the AND case, in a pipe with Poiseuille flow, we see generalists residing at the boundaries followed by specialists towards the middle. In the center where shear is lowest, cheaters are able to take over and destroy groups, leading to extinction. (c) For the OR case, in a pipe, we see coexistence between generalists and specialists at the boundary and specialists predominately in the center. The mutation rates used for the pipe are = 5 × 10 −7 s −1 for the AND case and = 5 × 10 −8 s −1 for the OR case. (b,d) In a Rankine vortex flow, we see generalists where shear is lowest, and specialists residing in an annulus where shear is at its maximum. The mutation rates used for the vortex are = 3 × 10 −7 s −1 for the AND case and = 5 × 10 −8 s −1 for the OR case. Secretion costs used are = 0.12 for the AND fitness, and = 0.17 for the OR fitness. The rest of the parameters are as given in Table 1. Average pipe population densities were obtained from averaging 5 runs in a domain of size 200 mm × 40 mm at a simulation time of = 2 × 10 6 s. Average vortex population densities were obtained from averaging 5 runs in a domain of size 60 mm × 60 mm also at a simulation time of = 2 × 10 6 s.

268
Specialization: • Large waste diffusion: Larger waste diffusion lowers self-competition and allows specialists to form denser groups to better utilize public goods secreted by neighbors. • Large public good benefit: A high benefit for public goods allows specialists to still be fit without secreting as many public goods. This also helps cheaters exploit producers. • Group structure: Specialists from groups when waste diffusion is larger than public good diffusion and when costs are not too low. When specialists do not form groups, they are easily taken over with cheaters, leading to either "chasing cheaters" (Figure 4-video 4) or extinction. When generalists form smaller, reproducing groups, they are able to escape take-over by specialists and out-compete specialists. Cheater coexistence: • Lack of group structure: Cheaters cannot exist on their own, but must "predate" on producers -generalists or specialists. When producers are fit, and do not form groups, they can grow quicker than cheaters fully taking over. This occurs when waste diffusion is large, and when secretion costs are low. • Small invasion fitness: When cheaters take over slower than producers reproduce, they are able to coexist. This happens when secretion costs are low, since the advantage of not secreting is lower (Figure 4-video 4). Extinction: • When cheaters take over: when the mutation rate and invasion fitness of cheaters large enough such that they take over groups faster than they fragment. This happens when public good benefit is large and waste diffusion is large, we see this in the right-middle regions of plots in • Enhanced specialization: Shearing flow can help specialist groups fragment quicker than generalist groups, and therefore transition a population to contain more specialists (Figure 3-video 2).  (Figure 5-video 1, Figure 5-video 2).

312
It is well known that spatial structure is key in the evolution of cooperation. By forming reproducing 313 groups, organisms are able to combat takeover by cheating mutants. In a similar vein, spatial 314 structuring can have a large effect on the evolution of specialization. 315 Previously we had found that mechanical factors involving diffusion and fluid shear can over- allow for functional specialization and multiple kinds of community structures. 321 We have studied a physically realistic advection-diffusion system describing microbial growth 322 and evolution. Through this first principles approach, we determined the effects of diffusion, public 323 good benefit and production cost, as well as flow patterns on the evolution of specialization. We 324 showed that specialization is stable when a transition to specialization results in smaller, faster 325 fragmenting groups or when it induces a discrete group structure from a spatially homogeneous 326 generalist state. 327 In nearly all studies of specialization, the basic mechanism is a trade off between different We study the Turing patterns of stable systems of generalists and specialists in both AND and OR fitness forms. 507 508 Generalist Turing analysis 509 For generalists, we linearize the following systeṁ 3 = ∇ 2 3 + 3 3 ( ) 1 = 1 ∇ 2 1 − 1 1 + 1 3 2 = 2 ∇ 2 2 − 2 2 + 2 3 3 = ∇ 2 3 − 3 + 3 .
We enforce 12 − − 2 < 0 for stability, otherwise the dense state will become unstable as the hill forms become saturated. Since the linear and constant coefficients are also negative, by Descartes rule of signs, in order to have any positive solution, we need 12 − 2 > 0. This implies we have up to two positive solutions, since there are two sign changes. In the OR case, we have * 3 given by solving the quadratic equation, 2 ( 12 − −2 ) 2 3 + 1 2 12 − 12 − 2 (2 + 12 ) 3 − 2 12 = 0 For this to have a positive solution, we require 2 12 − 12 − 2 (2 + 12 ) > 0. Plugging this into our linearized system, we get the eigenvalue equation, (− 2 + ) = Λ . Where, If we denote by ( ) = − 2 + , the characteristic equation for this system is now given as Now since (tr( )) 2 − tr( 2 ) = ( + + ) 2 − 2 + 2 + 2 > 0, tr( ) + tr( ) − 2tr( )tr( ) = 2( + 12 + ( + )) > 0, and we require (tr( )) 2 − tr( 2 ) > 0 from linear stability, the linear term is also positive. Therefore, by Descartes rule of signs, the requirement for Turing instability reduces to, det(− 2 + ) > 0, for some range of > 0. Inspecting this further, we see that det(− 2 + ) gives a cubic polynomial in 2 which goes to negative infinity as → ∞ and passes through det( ) < 0 at = 0. We therefore see that a Turing instability begins when the local maximum of Ψ( 2 ) = det(− 2 + ) vanishes. We therefore have that the critical wavenumber is given by maximizing Ψ( 2 ). Let = 2 , then The local maximum is then given by taking dΨ d = 0 and taking the positive root, where = 3det( ), = 2( , − + ,12 ), = det( ). Given the critical value * and a public good diffusion , we can find the critical waste diffusion that gives a Turing instability by setting Ψ( * , 3 ) = 0 and solving for 3 . This then gives a relation for 3 as a function of the other system parameters, when it exists. The theoretical region of Turing instability is compared against numerical results in Figure 2. The specialist case can be treated in a similar manner. Here we take 3 = 0 and because of symmetry, take 1 = 2 ≡ . We therefore again reduce our system to three equations and get the following system,̇ = ∇ 2 + ( ) 12 = ∇ 2 12 − 12 + ̇ 3 = ∇ 2 3 − 3 + 2 The factor of 2 in the waste chemical term comes from the fact that both populations 1 and 2 secrete waste while only one secretes the public good. We also have that the factor of 2 in front of the cost in the fitness functions now becomes a one. We therefore have the same results as in the generalist case with the substitutions, With these substitutions, the results from the generalists case applies also to the specialists. In general, the Turing domain will be different for the two states. These differing Turing regions and pattern types effect what type of state, generalist or specialist, will be more stable evolutionarily when also considering mutations. Effective group models 632 To get a better understanding of the results in Figure 3, we developed an effective model of ordinary differential equations describing the dynamics of each type of group. We assume that each type of species, generalists and specialists, form groups that grow and reproduce with different rates. The growth of the groups are given by a logistic equation, with a carrying capacity dependent on the domain. In order to fully describe the logistic growth, we also need to include transient states that are not stable but are continuously formed. The transient states therefore restrict the space in which the stable states are allowed to grow. We denote generalists groups by . Groups composed of generalists and specialists are transitioning to a specialist state and are denoted by . Mixed specialist groups made of microbes that secrete different public goods are represented by . Pure specialists groups composed of microbes secreting only one public good are denoted by . Pure specialist groups only occur in the OR case, and we exclude them in our effective model for the AND case. Finally, groups with cheaters that secrete no public good are denoted by . The evolutionary paths the groups can take are illustrated in Figure 1c. Effective model for AND type fitness 647 The group dynamics for the AND case are described by the following set of equations, where the constants correspond to the number of microbes in a group of type , and Σ = + + + is the total population. The rates give the reproduction rates for stable groups and the death rates for transient groups. The carrying capacity for the system is given by . Mutations cause groups to go down in number of secreted compounds (Figure 1c), and are given by terms proportional to the mutation rate . Back mutations do not stick in groups since they are paying more cost than their neighbors, and so we neglect transitions towards secreting more goods. The factor of 1/2 in front of mutation terms is due to mutations not always sticking in a group before it splits in two. The parameter values were obtained by fitting logistic growth curves to simulations without mutations to determine the natural growth of an isolated species. For this to be positive, the generalists group reproduction rate must satisfy,
We plot the results of our effective model as solid curves against simulation results (Figure 3) and get good overall agreement with the numerical simulations. At low costs ( ), microbes form stripes or become homogeneous, and the groups structure assumption of our effective model breaks down. Effective model for OR type fitness 686 In the OR case, pure specialist groups are now stable. Since we treat both chemicals symmetrically, and since the chemicals enter additively in the fitness function, pure specialist and mixed specialist groups are equivalent. However, to form a mixed specialist state, opposite mutations are needed within the lifetime of a generalist group. Since transition states are no longer unstable, once a mutation occurs in a generalist group, it quickly stabilizes to a specialist group. With these considerations, we build and solve an effective model only including generalists, pure specialists, and cheaters. The effective model is described by the following dynamical equations, Setting (Equation 9 -Equation 11) to zero to get steady states, we again get two non-trivial solutions. One with only specialist and cheating groups, and one where all three coexist. For the solution with only specialists and cheaters, we get = (2 − ) For the solution where all three coexist, we get for the generalists, * = (2 − )( − ) (2 − − 4 ) .