A Keller-Segel model for C elegans L1 aggregation
Leon Avery, Brian Ingalls, Catherine Dumur, Alexander Artyukhin
AA Keller-Segel model for
C elegans
L1 aggregation
Leon Avery , Brian Ingalls , Catherine Dumur , Alexander Artyukhin , Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario,Canada Department of Pathology, Virginia Commonwealth University, Richmond, VA, USA Chemistry Department, State University of New York, College of EnvironmentalScience and Forestry, Syracuse, NY, USA* [email protected] 8, 2020 1/39 a r X i v : . [ q - b i o . M N ] D ec bstract We describe a mathematical model for the aggregation of starved first-stage
Celegans larvae (L1s). We propose that starved L1s produce and respond chemotacticallyto two labile diffusible chemical signals, a short-range attractant and a longer rangerepellent. This model takes the mathematical form of three coupled partial differentialequations, one that describes the movement of the worms and one for each of thechemical signals. Numerical solution of these equations produced a pattern ofaggregates that resembled that of worm aggregates observed in experiments. We alsodescribe the identification of a sensory receptor gene, srh–2 , whose expression isinduced under conditions that promote L1 aggregation. Worms whose srh–2 gene hasbeen knocked out form irregularly shaped aggregates. Our model suggests thisphenotype may be explained by the mutant worms slowing their movement morequickly than the wild type.
Author summary
Among the most complex of animal behaviors are collective behaviors, in whichanimals interact with each other so as to produce large-scale organization. Starvedfirst-stage larvae of the nematode
Caenorhabditis elegans exhibit such a behavior: theycome together to form aggregates of several hundred worms. How and why they do thisare unknown. To address these questions, we developed a mathematical model ofstarved L1 aggregation. This model reproduced the main features of the behavior.
Introduction
Among the most complex behaviors exhibited by the nematode
C elegans are socialbehaviors such as mating [2] and aggregation. We recently described a new aggregationbehavior in starved
C elegans first-stage larvae (L1s) [1]. This new behavior raises twobroad questions whose answers we lack: (1) How do starved L1s aggregate? I.e., whatare the behavioral mechanisms by which they come together? (2) Why do starved L1saggregate? What selective advantage (if any) do these mechanisms or aggregation itselfprovide? To aid in answering these questions, we describe here a simple mathematicalDecember 8, 2020 2/39odel for L1 aggregation.L1 aggregation is not the only known
C elegans aggregation behavior, and ours isnot the first mathematical model of
C elegans aggregation. It has been known for manyyears that in the presence of food (bacteria), most true wild isolates of
C elegans aggregate, a behavior known as social feeding [3, 4]. Wild strains of
C elegans prefer lowconcentrations of oxygen. The usual
C elegans laboratory strain, N2, does not displaysocial feeding at normal atmospheric oxygen pressure because of a gain-of-functionmutation in the neuropeptide receptor gene npr–1 [3]. We and others have speculatedthat the consumption of oxygen in an aggregate of worms lowers oxygen concentrationand thereby attracts more worms [5], although this explanation is disputed [6]. Amathematical model of social feeding has recently been published [6].A third type of aggregation is mediated by indole-containing ascarosides [7]. L1s of daf–22 mutants, which are unable to make ascarosides [8, 9, 10]) aggregated similarly towild type [11], thus L1 aggregation is different from ascaroside-mediated aggregation.Observations of yet another type of aggregation have recently been published, togetherwith a model [12]. This form occurs in the long-term survival form of the worm—thedauer larva—and is probably mediated largely by a simple physical mechanism, surfacetension.The model we present here is simpler than previous
C elegans aggregation models inthe following sense: It does not describe aggregation behavior in completely realisticdetail. We attempt only to reproduce the essential aspects of the behavior. Accordingly,we simply assume the existence, which has been experimentally demonstrated [13, 14], oftaxis mechanisms that allow worms to move in the direction they want to go. Althoughtaxis mechanisms have been investigated for years and much is known about them [e.g.15–19], the model presented here is based on the idea that the end result of taxis(movement towards favored places) is sufficient to understand aggregation, and thatmechanistic details are not essential.A further simplification is to describe worms not as individuals, but via populationdensity: a continuous function of space and time, ρ ( t, x ). We also propose a simplemechanism for interactions among worms via diffusible chemical signals. The resultingmodel takes the form of a system of partial differential equations (PDEs), a variation onthe classic Keller-Segel [20] model, developed to explain the aggregation of cellular slimeDecember 8, 2020 3/39old amoebae. Our model has the advantage of high mathematical tractability, bothanalytical [21] and numerical, the latter of which is the focus of this paper. Results
Strategy
As described in the Introduction, our model is a deliberately simplified description ofbehavior that simply assumes the existence of taxis mechanisms that allow worms tomove in the direction they want to go. Further, worms are modeled not as individuals,but as a continuous function of space and time, population density ρ ( t, x ). Thissimplification allows us to describe worm movement with a mathematically tractablepartial differential equation (PDE) model similar to the well-studied Keller-Segel [20]model. Also, since worm density is a component of the model, it is straightforward todesign a model in which worm movement explicitly depends on density. Thedisadvantage is that individual worms are not accurately represented by a continuousdensity function. Moreover, we ignore the fact that worms are worm-shaped. I.e., aworm is long and thin, with a head at one end and a tail at the other. Worm geometryis central to some other published models of C elegans aggregation [5, 6, 12].We present two versions of the model: the precursor attractant-only model and thefinal attractant+repellent model. We begin with the simpler attractant-only model,which is closest in form to the original Keller-Segel model. This model was a partialsuccess—it reproduced certain aspects of L1 aggregation seen in experiments, whilefailing in others. The more complicated attractant+repellent model better reproducedL1 aggregation behavior.
Design of the PDEs
The attractant–only model for L1 aggregation consisted of two coupled PDEs, areaction-diffusion equation that describes the time evolution of the concentration of adiffusible chemical attractant, and a Fokker-Planck equation that describes themovement of the worms. The attractant PDE is:December 8, 2020 4/39 U = U t = − γU + D ∇ U + sρ (1)The worm PDE is: ˙ ρ = ρ t = ∇ · ( ρ ∇ ( V U ( U ) + V ρ ( ρ ) + σ log ρ )) (2)= ∇ ρ · ∇ V + ρ ∇ V + σ ∇ ρ (3)Functions and parameters appearing in the PDEs are ρ ( t, x ) > U ( t, x ) > σ > γ > D > s > V U ( U ) is a potential function that describes the worm’s response to attractant V ρ ( ρ ) is a potential function that describes the worms’ direct response to other worms V is a shorthand notation for V ( ρ, U ) = V U ( U ) + V ρ ( ρ )These equations are similar to those developed by Keller and Segel [20] to model theaggregation of Dictyostelium discoideum amoebae. Attractant PDE (1) is identical tothe reaction-diffusion equation with which they model acrasin. Worm PDE (2) is ageneralization of the equation they use to model the movement of amoebae, which, withspecific choices of the potential functions V U and V ρ , reduces to theirs.In designing this model for L1 aggregation, we sought to reproduce certain generalcharacteristics that were obvious in recordings of worm aggregation. First, the wormsaggregate. This suggests that they are somehow attracted to each other. Given what weknow about C elegans biology, it was an obvious guess that this attraction could bemediated by a diffusible chemical signal with limited range [14]. PDE (1) is essentiallyDecember 8, 2020 5/39he simplest physically plausible that meets these criteria.The design of the PDE describing the movement of the worms was more complicated.On the time scale of the experiments, there is neither birth nor death of new worms.This suggested that it should be possible to express the rate of change of worm densityas the divergence of some flow field. Since flux occurs by movement of worms, the netflux vector at any point is the density times the mean velocity of worms at that point.These considerations lead to a general equation of the form ρ t ( t, x ) = − ∇ · ( ρ ( t, x ) v ( ρ ( t, x ) , U ( t, x ) , ∇ ρ ( t, x ) , ∇ U ( t, x )))) (4)in which velocity v is a vector field depending on density and attractant concentrationand their gradients. We chose to assume that the velocity field is conservative, i.e. thatit can be represented as the gradient of some scalar potential field. There is nocompelling biological necessity for this assumption. We made it for two reasons: First itmakes the PDE system more tractable analytically. Second, in recordings of wormbehavior, we see that the worms eventually approach an equilibrium in which there islittle net flow of animals, and no cyclic flows are obvious.If the velocity field is conservative, then the velocity field v can be expressed as thenegative of the gradient of some scalar potential field V . We chose a potential functionthat is a sum of a signal-dependent potential V U and a density-dependent potential V ρ ,for convenience in separately engineering signal and density dependence. This led to thefinal form (2).For an attractant, V U must be a decreasing function of signal. In early simulationswith a linear V U numerical instability was a problem. Steep signal gradients frequentlyoccur in the course of simulation. With a linear V U , these led to large velocities, whichmeant that worm density at one location was rapidly affected by density at distantlocations. As a result it was impractical to satisfy the Courant–Friedrichs–Lewy (CFL)stability condition. We therefore sought a potential function whose dependence on U was convex. The Weber-Fechner Law, which is approximately true for many sensoryresponses [32], holds that response depends linearly on the logarithm of the stimulusintensity. (Keller and Segel [20] used a Weber-Fechner Law chemotactic responsefunction.) The Weber-Fechner Law has a serious defect: it implies that responseDecember 8, 2020 6/39ensitivity becomes infinite as concentration approaches zero. (I.e.,lim U → + d(log U )/d U = ∞ .) For these reasons we chose the following function. V U ( U ) := − β log( α + U ) (5)Parameter β determines the strength of attraction. At high signal concentrations (5)approximates a Weber-Fechner Law response. But at low signal concentrations, theresponse is linear. Parameter α determines where the transition from linear tologarithmic dependence takes place.We speculated that the circular shape of the aggregates [1] results from the wormspacking together as tightly as possible. To reproduce this effect in simulations, wedesigned a V ρ potential that would reflect worms taking up space. The ideal would havebeen a hard sphere potential V ρ ( ρ ) = ρ < ρ max ∞ if ρ ≥ ρ max (6)This potential function implies discontinuous time or spatial derivatives of density, andtherefore functions poorly with numerical methods for solving the PDE system. Wetherefore approximated the discontinuity with a hyperbolic tangent function. V ρ ( ρ ) = σ scale (cid:18) (cid:18) ρ − ρ max cushion (cid:19)(cid:19) (7)Four parameters determine the exact shape of V ρ : σ , ρ max , scale , and cushion . (Werefer to the latter two by the symbols used to represent them in software code, sincethey will play little role in the mathematics.) Two parameters, σ and scale , determinethe vertical scale. σ is the parameter that measures random worm movement (see (2)). V ρ rises from near 0 for small values of ρ to σ × scale for large ρ . ρ max is the density atwhich V ρ reaches half its maximum possible value. It is the point at which the V ρ curveis steepest, and therefore the closest approximation to ρ max of (6). cushion determineshow abrupt the rise of V ρ is. These functions are plotted in Fig 1. The Keller-Segelliterature describes other, less-flexible models in which organisms take up space [33],which we elected not to use.December 8, 2020 7/39 igure 1. Potential function plots. Potential functions that appear in the ρ PDE (2). Both potentials are made dimensionlessby dividing them by σ . Parameter values are as in Table 1. Parameter estimates
We required numerical estimates of parameters γ , D , and s that appear in (1) and σ of (2). In addition, we required ρ max , scale , cushion and α and β which determinethe shapes of the potential functions V ρ and V U .A C elegans
L1 is approximately a cylinder of diameter 15 µ m and length 240 µ m [34].December 8, 2020 8/39ince worms lie on their sides, a worm occupies approximate area 15 × ≈ µ m .We chose the inverse of this area, 28 000 cm − d as the parameter ρ max . ( d = 1 or 2 is thespatial dimension. The same number, 28 000, was used for one and two-dimensionalsimulations to facilitate comparison.) Parameters cushion and scale have no realbiological significance. Parameter cushion makes the ideal hard-sphere potential (6)continuous and differentiable, so that the PDEs can be solved numerically withdifferentiable functions. The value 2000 cm − d worked. The scale parameter need onlybe chosen large enough to constrain the maximum density—we chose 2.In small-scale simulations, we chose a mean density ¯ ρ = 9000 cm − d so thataggregates would occupy about 1 / / × − cm s − to 1 × − cm s − . (We assumed that the signals diffuse through theagar-solidified water under the worms. Worms also respond chemotactically to volatilechemicals diffusing through the air above them [14]—diffusion constants for suchvolatile signals would be much larger than for water-soluble signals.) We chose thediffusion constant of attractant, D a =1 × − cm s − , at the lower end of the range ofdiffusion constants in water. The aggregates that form have diameters of hundreds ofmicrometers. The mean distance a molecule of attractant diffuses before decaying is (cid:112) D a /γ a . We therefore chose γ a , the decay rate of attractant, to give it a range (cid:112) D a /γ a = 100 µ m. To fulfil its role in the model, repellent needs to have a longerrange, so we chose a large diffusion constant of D r =1 × − cm s − and a smallerdecay rate, giving it a range of 1 mm.Parameters s a and s r , the rates at which a worm secretes attractant and repellent,effectively set the units of concentration. We chose units of concentration such that s i and γ i (for i = a or r ) were numerically equal. (That is, if concentration is measured in“number of units of stuff”/cm d , we chose the units in which “stuff” is measured to be theamount secreted by one worm in one mean lifetime of the stuff, i.e. γ − i . This has theeffect that if γ i = (cid:104) number (cid:105) s − , then s i = (cid:104) number (cid:105) “stuff units” cm − d s − , with thenumber being the same in the two cases.) This ensures that concentrations U i andworm density ρ are in the same range numerically.Artyukhin et al found that the minimum worm density for aggregation is1500 cm − [1]. We identified this with the density threshold for instability. We choseDecember 8, 2020 9/39 a = α r = 1500 cm − , to make V U linear near the threshold, and to be obviouslyconvex near ρ max . We then chose β a = 2 σ to reproduce the 1500 cm − densitythreshold for instability in the attractant-only model. For the attractant+repellentmodel, we kept this value for β a and chose β r = − σ for the repellent. Adding repellentto the model increased the calculated instability threshold to 2357 cm − d .Parameter σ determines how rapidly the worms spread. Artyukhin et al [1] foundthat worms placed at the center of a 6 cm diameter petri plate spread to occupy muchof the area of the plate in 12 h, but at this time they still remain mostly concentratednear the center. We began by choosing values of σ that would approximately reproducethis distribution if the worms’ motions were purely diffusional, i.e., if their motions weregoverned by (8) ρ t = ˙ ρ = σ ∇ ρ, (8)with Neumann boundary condition d ρ d r (cid:12)(cid:12)(cid:12)(cid:12) r = R = 0 (9)Here R = 3 cm is the radius of the petri plate. Eqs. (8, 9) can be solved by separationof variables. Any solution ρ ( t, r, θ ) can be represented as a sum of exponentiallydecaying eigenfunctions of the Laplacian. The circularly symmetric eigenfunctions onthe disk with Neumann boundary condition (9) are ρ ( t, r ) = e − σk t J ( kr ) , (10)where J is a Bessel function of the first kind. The wavenumber k must be chosen tosatisfy boundary condition (9), i.e. k = j ,n /R , where j ,n is the n th nontrivial zero of J . The smallest wavenumber, corresponding to the circularly symmetric eigenfunctionthat decays most slowly, is thus k = j , /R ≈ . / ≈ .
28 cm − . We began bychoosing σ so that the corresponding time constant σ − j − , R was approximately 12 h.Of course, the motion of the worms is not purely diffusional. We therefore refined ourestimate of σ by numerical solution of the attractant+repellent system, from an initialcondition in which the worms began near the center of the petri plate. From theseDecember 8, 2020 10/39 able 1. Parameter values ¯ ρ mean worm density 9000 cm − d σ random worm movement 5 . × − cm s − ρ max midpoint of V ρ potential rise 28 000 cm − d cushion breadth of V ρ rise 2000 cm − d scale height of V ρ rise 2 β a strength of attraction 1 . × − cm s − α a attractant concentration scale 1500 cm − d γ a attractant decay rate 0 .
01 s − D a attractant diffusion constant 1 × − cm s − s a attractant secretion rate 0 .
01 cm − d s − β r strength of repulsion − . × − cm s − α r repellent concentration scale 1500 cm − d γ r repellent decay rate 0 .
001 s − D r repellent diffusion constant 1 × − cm s − s a repellent secretion rate 0 .
001 cm − d s − simulations we chose σ =5 . × − cm s − as producing results that resembledexperimental results.Table 1 summarizes parameter values. Simulation results, attractant-only model
Fig 2 shows the state of a numerical simulation of the attractant-only model after200 000 s (2 days and 7 hours). The initial condition was a uniform worm density of ¯ ρ =9000 cm − d , perturbed by normally distributed random noise of standard deviation 1%(i.e. 90 cm − d ). (The entire time courses can be seen in the videos options138a.mp4 and options139.mp4 in the Supporting Information.) Supporting Information Fig S1shows results at t = 200 000 s and t = 1 × s (116 days) of ten independent runs ofthe same simulation with different pseudorandom noise in the initial condition.This model successfully reproduced the experimental results in certain ways, butfailed in others. It was successful in that circular aggregates of maximum densityrapidly formed. The aggregates had sharp boundaries, and outside of aggregates wormdensity was low and uniform. This is most easily seen in the one-dimensional results(Fig 2A), but is also true in two dimensions. These are also characteristics of theexperimental results. (See Supporting Information video N2 5e5 washed.avi .)The model failed to reproduce the patterning of aggregates. In experiments(Supporting Information video
N2 5e5 washed.avi ), aggregation appears to reach anDecember 8, 2020 11/39 igure 2. Simulation of the attractant-only model
Panels
A, B show the results of simulations in one-dimensional space;
C, D show resultsin two-dimensional space.
A, C show density ρ ; B, D show attractant concentration U .The two numbers below each plot are the minimum and maximum values of the plottedfunction over the entire 1 cm × µ m in diameter. Most of the worms, in fact, end up in aggregates close tothis maximum size. These aggregates are also fairly uniformly spaced—the distancefrom one aggregate to its nearest neighbors varies little.In numerical solutions of the attractant-only model, however, aggregates had nomaximum size (aside from that imposed by the fixed finite number of worms), and theirspacing was not uniform. Furthermore, even after 200 000 s, they were not atequilibrium. This can be seen by computing the velocity v = (cid:107) ∇ V (cid:107) , which atequilibrium would be zero everywhere, but remained well above zero throughout thesimulation. More obviously, it is seen by continuing the solution past t = 200 000 s. FigS1 shows that aggregates increased in size and decreased in number between t =200 000 s and t = 1 × s. In fact, we believe the only true equilibria of theattractant-only model are those in which there is a single large aggregate containingalmost all the worms. This state was reached at t = 1 × s in one of the tensimulations in Fig S1.In fact, this observation is consistent with linear stability analysis of theattractant-only model (see Supplemental Information). PDE system (1, 2) shares withthe original Keller-Segel system the property of density-dependent instability. Thecondition for a sinusoidal variation of wavenumber k to be unstable is (30). Condition(30) has no nontrivial minimum wavenumber. That is, there is no natural maximumscale for the aggregates that form when the density exceeds threshold. It is true thatthe attractant has a natural range, (cid:112) D/γ —the distance an average molecule diffusesbefore it decays. However, worms attract each other, albeit weakly, even when they arefurther apart than this. There is thus no mechanism in the model (1, 2) that wouldprevent the merging of aggregates to unlimited size. This is true for any attractant-onlyKeller-Segel model.December 8, 2020 13/39 repellent is necessary
We could not reproduce the experimental observed scaling behavior in numericalexperiments with attractant-only models. By analogy with activator-inhibitor models ofbiological pattern formation, we suspected that the addition of a negative signal tooppose the attractant, a repellent, would solve the scale problem. Linear stabilityanalysis supports this intuition. (See Supporting information section Linear stabilityanalysis of the attractant+repellent model.)Intuitively, and by analogy with activator-inhibitor models, what one requires is ashort-range attractant and a long-range repellent. We therefore added to theattractant-only model a repellent with diffusion constant D r = 10 D a and decay rate γ r = 0 . γ a . The range of this repellent κ r = (cid:112) D r /γ r = 1 mm is ten times that ofattractant, and is approximately equal to the observed spacing between aggregates.Thus, in the attractant+repellent model, attractant PDE (1) is replaced with two PDEs,one (11) for attractant and the other (12) for repellent.˙ U a = − γ a U a + D a ∇ U a + s a ρ (11)˙ U r = − γ r U r + D r ∇ U r + s r ρ (12)As shown in Fig 3, addition of a repellent to the model produced the predicted effect.(Supporting Information videos options140a.mp4 and options141.mp4 show the fulltime courses for these simulations. Supporting Information Fig S2 shows tenindependent solutions of the two-dimensional system with different pseudorandom noiseat time 0.) Aggregates formed with characteristic size and spacing approximatelymatching those seen in experiments on worms. These solutions are close to equilibriumat t = 200 000 s, as seen by comparison with the results at t = 1 × s. (Compare FigsS2A and B.) There is even a hint of pattern formation, with the aggregates in anapproximate hexagonal array at t = 200 000 s. The hexagonal patterning is near perfectat t = 1 × s, with the exception of lattice defects, most easily recognized as slightlysmaller aggregates surrounded by five rather than six neighbors. If there is a problemwith the result, it is that the array is too perfect. More irregularity is seen inDecember 8, 2020 14/39 igure 3. Attractant+repellent simulationPanels
A, D show density ρ , B, E show attractant concentration U a , and C, F showrepellent concentration U r . The spatial units are centimeters. The two numbers beloweach plot are the minimum and maximum values of the plotted function over the entire1 cm × srh–2 and cooldown srh–2 encodes a G–protein coupled receptor expressed instarving L1s Attempting to understand molecular mechanisms of L1 aggregation, we measuredgene expression in starved L1s in the presence and absence of ethanol or acetate, one ofwhich is required for aggregation [1]. We identified an ethanol-induced gene, srh–2 ,whose expression increases at the time that starved L1s become capable of aggregation(Supporting information). srh–2 is predicted to encode a sensory receptor, i.e., a proteinexpressed on the surface of a sensory neuron, capable of detecting chemicals in theDecember 8, 2020 15/39nvironment. To find out if srh–2 plays a role in L1 aggregation, we knocked the geneout. (That is, we genetically engineered a mutant strain that lacks a functional srh–2 gene.) We then tested the srh–2 knockout worms for aggregation. As shown in Fig 4,these mutant worms still aggregate, but the aggregates are irregular in shape (Fig 4).
Figure 4. srh–2 knockout L1 aggregationStarved L1s of mutant worms lacking a functional srh–2 gene aggregate, but theaggregates they form are irregularly shaped (the animal crackers phenotype).
The Srh–2 phenotype may be modeled by rapid decay of wormmovement
Two observations suggested a partial explanation of the srh–2 phenotype. First, inthe movies of the attractant+repellent simulation, one sees formation of irregularlyshaped aggregates at early times. With time, these aggregates become circular. Second,in movies of aggregating L1s, there is a lot of rapid movement at early times, but as timegoes on, fewer worms are seen moving. This suggested that worm movement might slowwith time, perhaps because the worms run low on energy. (They are, after all, starving.)Together, these observations suggested an explanation for the Srh–2phenotype–—perhaps the movement of srh–2 knockout worms slows down faster. In themodel, such a movement slowdown would be reflected in the decrease of the parameters σ (representing random worm movement) and β a , β r (representing signal–responsivemovement) with time. We modeled slowdown with the attractant+repellent PDEs (2,11, 12), but with parameters σ, β a , β r time-dependent:December 8, 2020 16/39 = σ ( t ) = σ e − t/τ (13) β a = β a ( t ) = 2 σ ( t ) (14) β r = β r ( t ) = − σ ( t ) (15)(Note that this is not the same as simply stretching the time axis of theattractant+repellent model, since the time–scales of Eqs (11, 12) remain unchanged.)The t = 0 values σ (0) , β a (0) , β r (0) were the same as those of σ, β a , β r in Table 1. Fig 5shows results at t = 200 000 s of simulations of the slowdown model with four differentvalues of τ . When τ is very small (e.g. 30 min, Fig 5A) aggregation is arrested beforedense aggregates form. Larger values of τ permit the formation and persistence ofirregular aggregates. Full-scale simulations
Our experimental studies of aggregation usually begin by placing a large number ofworms in the center of a 6 cm petri plate [11]. (See Supporting video
N2 5e5 washed.avi which records the behavior over 12 h of 500 000 worms that wereplaced on the center of a plate at time 0.) These experiments begin with the dispersalof the worms, so that the density is high near the center of the plate and lower towardsthe edges. To more closely mimic such experiments, we solved the attractant+repellentmodel on a 6 cm × t = 200 000 s. (Supporting Information video options157.mp4 shows theentire time course.)We do not believe that the attractant+repellent model accurately represents thephysics and biology of worm motions in the region near the center of the plate. Forinstance, in the preparation of eggs from which the L1s used for the experiment hatch,some non-living debris is inevitably generated. This debris, which is transferred to thecenter of the plate along with the worms, may influence behavior.Outside this central region, the behavior of the full-scale simulation resembled theDecember 8, 2020 17/39 igure 5. Attractant+repellent simulation with slowdownWorm density ρ ( t, x ) of a slowdown model simulation at t = 200 000 s for four differentvalues of τ .behavior of worms on a 6 cm diameter petri plate. Discussion
Our final model, the attractant+repellent model, appeared to reproduce the mainfeatures of L1 aggregation. Further work is required to quantify how well the modelreproduces experimental results [21]. This model is minimal, we believe, in the senseDecember 8, 2020 18/39 igure 6.
Full-scale attractant+repellent simulationSimulation of the attractant+repellent model on a 6 cm × × srh–2 is expressed under conditions where L1 aggregation takes place. Mutantworms whose srh–2 gene has been knocked out aggregate, but their aggregates areirregularly shaped, unlike the uniformly circular aggregates of wild-type worms. Ourmodeling suggests this phenotype could be explained by a faster-than-normal decreasein worm movement in the mutant. This observation is potentially testable by trackingthe movement of individual fluorescently labeled worms. Materials and methods
Numerical solution of partial differential equations (PDEs)
We simulated
C elegans
L1 aggregation by solving PDEs (1, 2) or (11, 12, 2)numerically in one or two spatial dimensions. In models with repellent and attractant,there were three PDEs, one for ρ (2) and one each for U a (11), and U r (12) . Thedomain for one-dimensional simulations was a simple interval Ω = [0 , w ]. Domains fortwo-dimensional simulations were rectangular Ω = [0 , w ] × [0 , h ]. Width w and height h varied according to the problem. To avoid distortion of the behavior by boundaryeffects, all simulations were carried out with periodic boundary conditions. (For anexplanation and examples of these boundary effects, see Avery [21]).Continuous fields ρ , U a and U r were approximated by a grid of points equally spacedin each dimension. The spatial derivatives in the PDEs were replaced with linearDecember 8, 2020 20/39ombinations of the function values ρ , U a , U r , and V ( ρ, U a , U r ) (5, 7) to approximatethe time derivative of each field at each point to fourth order. Simulation of theattractant+repellent model at a resolution of 384 cm − on a 6 cm × × (6 × = 15 925 248 degrees of freedom.We implemented the solution of this system of ODEs (ordinary differentialequations) in PETSc (the “Portable Extensible Toolkit for Scientific computation”)[22, 23, 24]. Among the tools included in PETSc is the TS (time-stepper) package, alibrary of ODE/DAE (differential algebraic equation) solvers [25]. All solutions shownwere produced with the PETSc Rosenbrock-W time stepper ra34pw2 [26], an implicitthird-order method. We used PETSc’s basic adaptive step size mechanism. Thismethod uses error estimates from the embedded stepper to adjust step size so as tomaintain error below predetermined absolute and relative tolerances. In addition, weimposed a step size limit inspired by the Courant-Friedrichs-Lewy (CFL) condition. Ateach step we calculated the mean worm velocity v = − ∇ V at each point. We limitedstep size to min( | ∆ x/v x | , | ∆ y/v y | ). (Here ∆ x and ∆ y are the point spacing in the x and y directions, and the minimum is taken over both dimensions and all spatial points.)Linear equations were solved with the MUMPS parallel direct solver [27, 28] forone-dimensional problems and with PETSc’s built-in gmres (generalized minimalresidual) Krylov solver for two-dimensional problems. Full-scale simulations
Full-scale simulations were carried out on a 6 cm × ρ = 2000 cm − .(This mean density is much lower than the mean density ¯ ρ = 9000 cm − used forsmall-scale simulations with uniform density, because in the full-scale simulations wormsare concentrated near the center of the domain, where the density is higher than themean.) These 72 000 worms were made up of 3600 distributed uniformly on the plate toavoid zero or negative densities, which would result in (2) becoming undefined, plus68 400 worms placed in a 2 cm diameter circle at the center of the square, with thedensity in the square as if a 2 cm diameter sphere had been placed in the center and theworms in it fell vertically onto the surface.December 8, 2020 21/39 (0 , x, y ) = b ρ + a ρ (cid:115) max (cid:18) , − ( x − + ( y − R (cid:19) (16) a ρ = 3(¯ ρ − b ρ ) w πR (17) b ρ = 100 cm − (18) R = 1 cm (19)¯ ρ = 2000 cm − (20) w = 6 cm (21)To this initial condition was added normally distributed pseudorandom noise withstandard deviation 0 . ρ (0 , x i , y j ) at each grid point ( x i , y j ). In addition, to simulatethe continuous production of noise that occurs in real experiments, pseudorandom noisewas injected in the course of the simulation. Noise generation was modeled as anindependent geometric Brownian motion attached to each grid point. Noise was injectedat times t n = 10 n/ s for n from 0 to 10, i.e., at 1 s, 3 .
16 s, 10 s, 31 . ρ ( t n , x i , y j ) wasmultiplied by a pseudorandom number exp { P n ( x i , y j ) } where P n ( x i , y j ) areindependent normal pseudorandom variates with variance 10 − ∆ t , ∆ t being the amountof time passed since the last noise injection. After noise injection, ρ ( t n , x i , y j ) werenormalized so that the total number of worms remained unchanged. The timing of noiseinjection was a pragmatic compromise. Computational efficiency precludes injectingnoise continuously, since time steps had to be small immediately after noiseinjection—high spatial frequency noise resulted in rapid worm movement. Wormsmoved rapidly near the beginning of the simulation and more slowly near the end, asthey approach a stable equilibrium. We thus chose a schedule in which the frequency ofnoise injection decreased steadily with time.December 8, 2020 22/39 ffect of ethanol and acetate on the transcriptome of starvedL1s To obtain L1s for transcriptomic profiling, we grew N2 worms in liquid culture. Weinoculated 250 mL S-complete in a 2 L flask with 7 × synchronized L1 larvaeobtained from a small-scale liquid culture and added 10 mL 50% E. coli
K-12 stocksuspension. Worms were grown at 22 ° C, 220 rpm (for details see Artyukhin et al [1]).We monitored the worm culture during the next 2 . E. coli food as itbecame depleted. Bleaching of gravid adult worms (8 mL water + 2 mL bleach +0 . × eggs. After 3washes with M9 buffer, eggs were resuspended in 20 mL M9 and allowed to hatch. After2 days of starvation (51 h to 53 h after bleaching, 20 ° C), we collected L1 larvae, washedthem 6 times with M9, 10 mL each time, and resuspended in 10 mL M9 after the finalwash. 300 µ L of the resulting L1 suspension (corresponding to ca. 20 µ L L1 pellet) wereadded to each of the following solutions in 15-ml plastic tubes: (1) 3 mL M9; (2) 3 mLM9 + 3 µ L ethanol; (3) 3 mL M9 + 51 µ L 1 M potassium acetate. Tubes were put on arocker at room temperature (ca. 21 ° C). After 1 . . µ L Trizol to each sample,and froze in liquid nitrogen. The above procedure was repeated two more times withworms grown on different days to obtain biological triplicates for each condition(control, ethanol, potassium acetate). All samples were stored at -80 ° C before analysis.
RNA extraction
Total RNA was extracted and the quality evaluated using a sample processingmethod we described previously [29]. Total RNA was extracted from 400
C elegans worms using the MagMAX ™ -96 for Microarrays Total RNA Isolation Kit (Invitrogen ™ Life Technologies, Carlsbad, CA), in an automated fashion using the magnetic particleprocessors MagMAX ™ Express. RNA purity was judged by spectrophotometry at260 nm, 270 nm and 280 nm. RNA integrity as well as cDNA and cRNA synthesisproducts were assessed by running 1 µ L of every sample in RNA 6000 Nano LabChips ™ December 8, 2020 23/39n the 2100 Bioanalyzer (Agilent Technologies, Foster City, CA).
Gene expression microarray analyses
The Affymetrix ™ protocol utilized for our microarray analyses has been previouslydescribed [29], and was used with the following modifications. Starting from 500 ng oftotal RNA, we performed a single-strand cDNA synthesis primed with a T7–(dT24)oligonucleotide. Second strand cDNA synthesis was performed with the E. coli
DNAPolymerase I, and biotinylation of the cRNA was achieved by in vitro transcription(IVT) reaction using the using the GeneChip ™
3’ IVT Express Kit (Affymetrix, SantaClara, CA). After a 37 ° C incubation for 16 h, the labeled cRNA was purified using thecRNA cleanup reagents from the GeneChip ™ Sample Cleanup Module. As per theAffymetrix ™ protocol, 10 µ g of fragmented cRNA were hybridized on the GeneChip ™ Celegans
Genome array (Affymetrix Inc., Santa Clara, CA) for 16 h at 60 rpm in a 45 ° Chybridization oven. The GeneChip ™ C elegans
Genome array provides a comprehensivecoverage of the transcribed
C elegans genome by analyzing the expression level of over22 500 well-characterized transcripts. The arrays were washed and stained withstreptavidin phycoerythrin (SAPE; Molecular Probes, Eugene, OR) in the Affymetrix ™ fluidics workstation. Every chip was scanned at a high resolution, on the Affymetrix ™ GeneChip Scanner 3000 7G according to the GeneChip ™ Expression Analysis TechnicalManual procedures (Affymetrix). After scanning, the raw intensities for every probewere stored in electronic files (in .DAT and .CEL formats) by the GeneChip ™ OperatingSoftware v1.4 (GCOS) (Affymetrix). Overall quality of each array was assessed bymonitoring the 3 (cid:48) /5 (cid:48) ratios for the housekeeping gene, glyceraldehyde 3-phosphatedehydrogenase (Gapdh), and the percentage of “Present” genes (%P). Arrays exhibitingGapdh 3 (cid:48) /5 (cid:48) < . >
40% were considered good quality arrays.
Array normalization
For the microarray data analysis, background correction, normalization, andestimation of probe set expression summaries was performed using the log-scale robustmulti-array analysis (RMA) method [30]. Hierarchical cluster analyses were performedwith the BRB-ArrayTools v3.1.0 (Biometric Research Branch, National CancerDecember 8, 2020 24/39nstitute), an Excel add-in that collates microarray data with sample annotations. Inorder to identify differentially expressed genes between the different classes, weperformed t-tests for each probe set from biological replicates in each class. Statisticalsignificance for multivariate analysis to assess probe set specific false discovery rates(FDR) was performed by estimating the q-values, using the Bioconductor q-valuepackage [31]. srh–2 knockout
Deletion mutants of srh–2 were generated by CRISPR and obtained from Knudra.Three deletion strains were generated (COP-1274, 1275, 1276), all of the same genotype: unc–119(ed3) III; srh–2(knu317::unc–119(+)) V . References
1. Artyukhin AB, Yim JJ, Cheong MC, Avery L. Starvation-induced collectivebehavior in C. elegans. Scientific reports. 2015;5:10647. doi:10.1038/srep10647.2. Barr MM, Garcia LR. Male mating behavior.; 2006.3. de Bono M, Bargmann CI. Natural variation in a neuropeptide Y receptorhomolog modifies social behavior and food response in C. elegans. Cell.1998;94(5):679–689.4. de Bono M, Tobin DM, Davis MW, Avery L, Bargmann CI. Social feeding inCaenorhabditis elegans is induced by neurons that detect aversive stimuli. Nature.2002;419(6910):899–903.5. Rogers C, Persson A, Cheung B, de Bono M. Behavioral motifs and neuralpathways coordinating O2 responses and aggregation in C. elegans. Curr Biol.2006;16(7):649–659. doi:S0960-9822(06)01316-9 [pii] 10.1016/j.cub.2006.03.023.6. Ding SS, Schumacher LJ, Javer AE, Endres RG, Brown AE. Shared behavioralmechanisms underlie C. elegans aggregation and swarming. eLife. 2019;8.doi:10.7554/eLife.43318.December 8, 2020 25/39. Von Reuss SH, Bose N, Srinivasan J, Yim JJ, Judkins JC, Sternberg PW, et al.Comparative metabolomics reveals biogenesis of ascarosides, a modular library ofsmall-molecule signals in C. elegans. Journal of the American Chemical Society.2012;134(3):1817–1824. doi:10.1021/ja210202y.8. Butcher RA, Ragains JR, Li W, Ruvkun G, Clardy J, Mak HY. Biosynthesis ofthe Caenorhabditis elegans dauer pheromone. Proc Natl Acad Sci U S A.2009;106(6):1875–1879. doi:0810338106 [pii] 10.1073/pnas.0810338106.9. Golden JW, Riddle DL. A gene affecting production of the Caenorhabditiselegans dauer-inducing pheromone. Mol Gen Genet. 1985;198(3):534–536.10. Joo HJ, Yim YH, Jeong PY, Jin YX, Lee JE, Kim H, et al. Caenorhabditiselegans utilizes dauer pheromone biosynthesis to dispose of toxic peroxisomalfatty acids for cellular homoeostasis. The Biochemical journal. 2009;422(1):61–71.doi:10.1042/BJ20090513.11. Artyukhin AB, Yim JJ, Cheong Cheong M, Avery L. Starvation-inducedcollective behavior in C. elegans. Scientific Reports. 2015;5.doi:10.1038/srep10647.12. Sugi T, Ito H, Nishimura M, Nagai KH. C. elegans collectively forms dynamicalnetworks. Nature Communications. 2019;10(1):683.doi:10.1038/s41467-019-08537-y.13. Bargmann CI. Chemosensation in C. elegans. WormBook. 2006; p. 1–29.doi:10.1895/wormbook.1.123.1.14. Bargmann CI, Mori I. Chemotaxis and Thermotaxis. In: Riddle DL, BlumenthalT, Meyer BJ, Preiss JR, editors. C. elegans II. New York: Cold Spring HarborPress; 1997. p. 717–737.15. Mori I. Genetics of chemotaxis and thermotaxis in the nematode Caenorhabditiselegans. Annu Rev Genet. 1999;33:399–422.16. Dunn NA, Lockery SR, Pierce-Shimomura JT, Conery JS. A Neural NetworkModel of Chemotaxis Predicts Functions of Synaptic Connections in theDecember 8, 2020 26/39ematode Caenorhabditis elegans. Journal of Computational Neuroscience.2004;17(2):137–147. doi:10.1023/B:JCNS.0000037679.42570.d5.17. Roberts WM, Augustine SB, Lawton KJ, Lindsay TH, Thiele TR, Izquierdo EJ,et al. A stochastic neuronal model predicts random search behaviors at multiplespatial scales in C. elegans. eLife. 2016;5. doi:10.7554/eLife.12572.18. Pierce-Shimomura JT, Morse TM, Lockery SR. The fundamental role ofpirouettes in Caenorhabditis elegans chemotaxis. The Journal of neuroscience :the official journal of the Society for Neuroscience. 1999;19(21):9557–69.doi:10.1523/JNEUROSCI.19-21-09557.1999.19. Ferr´ee TC, Lockery SR. Computational rules for chemotaxis in the nematode C.elegans. Journal of computational neuroscience. 2020;6(3):263–77.doi:10.1023/a:1008857906763.20. Keller EF, Segel LA. Initiation of slime mold aggregation viewed as an instability.Journal of Theoretical Biology. 1970;26(3):399–415.doi:10.1016/0022-5193(70)90092-5.21. Avery, Leon. Mathematical Modeling of C elegans L1 aggregation [Ph.D. thesis].University of Waterloo; 2020. Available from: http://hdl.handle.net/10012/15480 .22. Balay S, Abhyankar S, Adams MF, Brown J, Gropp P, Buschelman K, et al..PETSc Web Page; 2017.23. Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, et al. { PETS } c Users Manual. Argonne National Laboratory; 2019. ANL-95/11 -Revision 3.11. Available from: .24. Balay S, Gropp WD, McInnes LC, Smith BF. Efficient Management ofParallelism in Object Oriented Numerical Software Libraries. In: Arge E, BruasetAM, Langtangen HP, editors. Modern Software Tools in Scientific Computing.Birkh¨auser Press; 1997. p. 163–202.December 8, 2020 27/395. Abhyankar S, Brown J, Constantinescu EM, Ghosh D, Smith BF, Zhang H.PETSc/TS: A Modern Scalable ODE/DAE Solver Library. arXiv preprintarXiv:180601437. 2018;.26. Rang J, Angermann L. New Rosenbrock W-methods of order 3 for partialdifferential algebraic equations of index 1. BIT Numerical Mathematics.2005;45(4):761–787.27. Amestoy PR, Guermouche A, L’Excellent JY, Pralet S. Hybrid scheduling for theparallel solution of linear systems. Parallel Computing. 2006;32(2):136–156.28. Amestoy PR, Duff IS, Koster J, L’Excellent JY. A Fully AsynchronousMultifrontal Solver Using Distributed Dynamic Scheduling. SIAM Journal onMatrix Analysis and Applications. 2001;23(1):15–41.29. Dumur CI, Nasim S, Best AM, Archer KJ, Ladd AC, Mas VR, et al. Evaluationof Quality-Control Criteria for Microarray Gene Expression Analysis. ClinicalChemistry. 2004;50(11):1994–2002. doi:10.1373/clinchem.2004.033225.30. Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP. Summaries ofAffymetrix GeneChip probe level data. Nucleic acids research. 2003;31(4):e15.doi:10.1093/nar/gng015.31. Storey JD. A direct approach to false discovery rates. Journal of the RoyalStatistical Society Series B: Statistical Methodology. 2002;64(3):479–498.doi:10.1111/1467-9868.00346.32. Fechner GT. Elemente der Psychophysik. Leipzig: Druck and Verlag vonBreitkopf and H¨artel; 1860.33. Hillen T, Painter KJ. A user’s guide to PDE models for chemotaxis. Journal ofMathematical Biology. 2009;58(1-2):183–217. doi:10.1007/s00285-008-0201-3.34. Avery L, Shtonda BB. Food transport in the C. elegans pharynx. J Exp Biol.2003;206(Pt 14):2441–2457.35. Luca M, Chavez-Ross A, Edelstein-Keshet L, Mogilner A. Chemotactic signaling,microglia, and Alzheimer’s disease senile plaques: Is there a connection? BulletinDecember 8, 2020 28/39f Mathematical Biology. 2003;65(4):693–730.doi:10.1016/S0092-8240(03)00030-2.36. Wolfram Research I. Mathematica; 2019. Supporting information
Linear stability analysis of the attractant-only model
The attractant–only model (1, 2) shares with the original Keller-Segel system [20]the property of density-dependent instability. There is a uniform equilibrium, ρ eq ( t, x ) = ¯ ρ (22) U eq ( t, x ) = ¯ U := s a γ a ¯ ρ (23)Substituting these functions into (1, 2) shows that d ρ eq /d t = d U eq /d t = 0. Consider apopulation near this equilibrium, and let ρ ( t, x ) = ¯ ρ + δρ ( t, x ), U ( t, x ) = ¯ U + δU ( t, x ).Because (¯ ρ, ¯ U ) is an equilibrium, ρ t ( t, x ) = δρ t ( t, x ) and U t ( t, x ) = δU t ( t, x ).Substituting into (1, 2), δρ t ( t, x ) = ∇ · (cid:0) δρV (cid:48) U ( U ) ∇ ¯ U + ¯ ρV (cid:48) U ( U ) ∇ δU (cid:1) (24)+ ∇ · (cid:0) δρV (cid:48) ρ ( ρ ) ∇ ¯ ρ + ¯ ρV (cid:48) ρ ( ρ ) ∇ δρ (cid:1) + ∇ · ( δρV (cid:48) U ( U ) ∇ δU )+ σ ∇ ¯ ρ + σ ∇ δρδU t ( t, x ) = − γδU + D ∇ δU + sδρ (25) ∇ ¯ U = ∇ ¯ ρ = ∇ ¯ ρ = 0. Writing V (cid:48) U ( U ) = V (cid:48) U ( ¯ U ) + (cid:79) ( δU ), V (cid:48) ρ ( ρ ) = V (cid:48) ρ (¯ ρ ) + (cid:79) ( δρ ) andignoring second order terms we have, to first order, the linear vector PDE, ∂∂t δρδU = ( σ + V (cid:48) ρ (¯ ρ )) ∇ V (cid:48) U ( ¯ U )¯ ρ ∇ s − γ + D ∇ δρδU (26)December 8, 2020 29/39his is easily solved by separation of variables. It will be convenient to define σ (cid:48) = σ + V (cid:48) ρ (¯ ρ ). Since V ρ is, by design, an increasing function, V (cid:48) ρ > σ (cid:48) > σ > V (cid:48) ρ (¯ ρ ) ≈ σ (cid:48) ≈ σ . Now, eigenfunctions of (26) are ofthe form δρ ( t, x ) δU ( t, x ) = ρ k e λ k t e i k · x U k e λ k t e i k · x (27) k is a wavenumber vector, i.e., a vector of frequency in each spatial direction.Substituting into (26) produces the 2 × λ k ρ k U k = − σ (cid:48) k − V (cid:48) U ( ¯ U )¯ ρk s − γ − Dk ρ k U k (28)where k := (cid:107) k (cid:107) . Solutions of the linearized system (28) are linear combinations offunctions (27) where ( ρ k , U k ) (cid:124) is an eigenvector of the matrix in (28). If, for every k ,Re( λ k ) < k such that thematrix has an eigenvalue with positive real part, then the uniform equilibrium isunstable. The sum of the two eigenvalues, the trace of the matrix, is negative, − σ (cid:48) k − γ − Dk <
0, so the only possible way to have an eigenvalue with positive realpart is if both eigenvalues are real, one positive and one negative. Thus, first-orderinstability is expected if and only if the determinant is negative, (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − σ (cid:48) k − V (cid:48) U ( ¯ U )¯ ρk s − γ − Dk (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) < Dσ (cid:48) k + (cid:0) γσ (cid:48) + V (cid:48) U ( ¯ U ) s ¯ ρ (cid:1) k < γσ (cid:48) + V (cid:48) U ( ¯ U ) s ¯ ρ <
0. Since, as mentioned above, V U is a decreasing function of attractant concentration, V (cid:48) U ( ¯ U ) <
0. If the condition γσ (cid:48) + V (cid:48) U ( ¯ U ) s ¯ ρ < k , the negative k term willdominate the positive k term and the determinant will be negative. Thus, first-orderinstability occurs if and only if ¯ ρ > − γσ (cid:48) sV (cid:48) U ( ¯ U ) . Since ¯ U := s ¯ ρ/γ depends on ¯ ρ , it is not aDecember 8, 2020 30/39oregone conclusion that instability is possible. Neglecting V (cid:48) ρ and with V U as in (5), theinstability condition reduces to ¯ ρ > αγσs ( β − σ ) > β > σ . By design, the instability condition is ¯ ρ > − d with the parameter values in Table 1. Linear stability analysis of the attractant+repellent model
For simplicity, we assume V ρ = 0 in the following analysis. With the parametervalues in Table 1, this is very close to true in the vicinity of the instability threshold.Linearization of the PDE system (2, 11, 12) around a uniform equilibrium at ρ eq ( t, x ) = ¯ ρ , U i, eq ( t, x ) = ¯ U i = s i ¯ ρ/γ i (for i ∈ a, r ) produces the following linear PDEsystem dd t δρδU a δU r = σ ∇ ¯ ρV (cid:48) U a ( ¯ U a ) ∇ ¯ ρV (cid:48) U r ( ¯ U r ) ∇ s a − γ a + D a ∇ s r − γ r + D r ∇ δρδU a δU r (32)The ansatz δρδU a δU r = ρ k U a k U r k e λ k t e i k · x (33)yields the eigenvalue problemDecember 8, 2020 31/39 k ρ k U a k U r k = − σk − ¯ ρV (cid:48) U a ( ¯ U a ) k − ¯ ρV (cid:48) U r ( ¯ U r ) k s a − γ a − D a k s r − γ r − D r k ρ k U a k U r k (34)= − σk − ¯ ρk ( V (cid:48) ( ¯U )) (cid:124) ¯ ρ s − Γ − D k ρ k U k (35)= − N (¯ ρ, k ) ρ k U k (36)In (35), the matrix is in a block form that can be extended easily to any number ofsignals, with s = s a s r (37) U k = U a k U r k (38) V (cid:48) ( ¯U ) = V (cid:48) U a ( ¯ U a ) V (cid:48) U r ( ¯ U r ) (39) Γ = γ a γ r (40) D = D a D r (41)In (36), N (¯ ρ, k ) is defined as the negative of the matrix in (35). (We define N as thenegative to avoid an inconvenient factor of ( − n signals in the determinant we areabout to calculate.) The uniform equilibrium is unstable at mean density ¯ ρ if for some k , N (¯ ρ, k ) has an eigenvalue with negative real part. If (cid:12)(cid:12) N (¯ ρ, k ) (cid:12)(cid:12) <
0, the equilibriumis certainly unstable. This leads to the following criterion for instabilityDecember 8, 2020 32/39 σ > ¯ ρ ( V (cid:48) ( ¯U )) (cid:124) ( Γ + D k ) − s (42)= ¯ ρ (cid:88) i ∈{ a,r } V (cid:48) U i ( ¯ U i ) s i γ i + D i k (43)= ¯ ρ (cid:18) V (cid:48) U a ( ¯ U a ) s a γ a + D a k + V (cid:48) U r ( ¯ U r ) s r γ r + D r k (cid:19) (44)It is possible to choose parameter values so that this criterion predicts instabilitywith a nontrivial minimum wavenumber (and therefore finite maximum scale). Howdoes this work? Remember that V (cid:48) U r > V (cid:48) U a < γ r < γ a and D r > D a , because the repellent is a longer-range signal than the attractant. Thus, atlow k , if the relative magnitudes of V (cid:48) U r s r and V (cid:48) U a s a are appropriately adjusted (byevolution or the modeler), the sum in (44) is positive and the uniform equilibrium isstable to perturbations of small wavenumber = large scale. As k rises the D r k factorin the denominator of the repellent term makes the positive term small compared to thenegative attractant term. The sum in parentheses can become negative, and if ¯ ρ is largeenough, the right-hand-side drops below − σ , and instability to perturbations ofintermediate wavenumber = medium scale results. For large k the right-hand-sideapproaches zero because of the Dk factors in both denominators. The uniformequilibrium is thus stable to perturbations of large wavenumber = small scale. It istherefore possible for an attractant+repellent Keller-Segel model to have a finite naturalscale.Based on calculations of this sort we chose β r = − β a = − σ = − . × − cm s − .Attractant parameters remained as in Table 1. The addition of a repellent increases thethreshold for instability, but the uniform equilibrium is still predicted to be unstable at¯ ρ = 9000 cm − d. Convergence tests
To test whether the numerical solution of the worm system PDEs (2, 11, 12)approximates the correct solution, we compared numerical solutions of theattractant+repellent system to an analytical solution of the linearized PDE system (32).December 8, 2020 33/39n one dimension, for x ∈ Ω = [0 , δρ ( t, x ) := a ρ e λt sin( φ + 2 πk x ) (45) ρ lin ( t, x ) = ¯ ρ + δρ ( t, x ) (46) U a, lin ( t, x ) = s a γ a ¯ ρ + U a,k δρ ( t, x ) (47) U r, lin ( t, x ) = s r γ r ¯ ρ + U r,k δρ ( t, x ) (48)In two dimensions, for ( x, y ) ∈ Ω = [0 , √ / × [0 , / δρ ( t, x, y ) := a ρ e λt (cid:18) (cid:0) cos(2 πk y )+ (49)cos(2 π ( k x, x − k y, y ))+cos(2 π ( k x, x + k y, y )) (cid:1)(cid:19) ρ lin ( t, x, y ) = ¯ ρ + δρ ( t, x, y ) (50) U a, lin ( t, x, y ) = s a γ a ¯ ρ + U a,k δρ ( t, x, y ) (51) U r, lin ( t, x, y ) = s r γ r ¯ ρ + U r,k δρ ( t, x, y ) (52)Here φ, a ρ ∈ R and k ∈ Z are parameters that can be freely chosen. We chose a ρ = 1and φ = π/
2. We chose k = 4 to produce a substantial positive growth rate. For thetwo-dimensional case, we chose k x, = k √ / k y, = k / λ is the positive eigenvalue of matrix − N (¯ ρ, k ) (36), with correspondingeigenvector (1 , U ak , U rk ) (cid:124) (as in (33), but normalized so that ρ k = 1). Numericalvalues λ ≈ .
000 955 s − , U ak ≈ . U rk ≈ .
121 were estimated to 15-digitprecision by numerical diagonalization of the computed matrix − N (¯ ρ, k ).Functions ( ρ lin ( t, x ) , U a, lin ( t, x ) , U r, lin ( t, x )) (cid:124) are of course not an exact solution ofthe full nonlinear PDEs (2, 11, 12). To produce a closely related system with this exactsolution for convergence testing, we modified the ρ PDE (2) by addition of a sourceterm S ( t, x ).December 8, 2020 34/39 ρ = ∇ · ( ρ ∇ ( V U a ( U a ) + V U r ( U r ) + V ρ ( ρ ) + σ log ρ )) + S ( t, x ) (53) S ( t, x ) = ∂ρ lin ( t, x ) ∂t − ∇ · (cid:0) ρ lin ( t, x ) (54) ∇ (cid:0) V U a ( U a, lin ( t, x )) + V U r ( U r, lin ( t, x ))+ V ρ ( ρ lin ( t, x ))+ σ log ρ lin ( t, x ) (cid:1)(cid:1) It was unnecessary to modify the U a and U r PDEs since they are linear. Sourcefunction (54) was computed symbolically from linear solutions (46-48) or (50-52) andconverted to sympy expressions with Mathematica [36]. sympy is a computer algebrapackage for the programming language python . Software
Software developed for this work is available from https://github.com/leonavery/KSFD . Supplemental figuresSupplemental data
The following data file is provided:
C.elegans microarray results 020615 CID.xlsx
Microarray expression profilingresults.
Videos
The following video files are provided atAvery, Leon (2020), “Avery L1agg2”, Mendeley Data, V1, doi: 10.17632/r5v772ftcs.1http://dx.doi.org/10.17632/r5v772ftcs.1
N2 5e5 washed.avi
This video shows the time course of aggregation after 500 000starved L1s were pipetted onto the center of a petri plate. The recording coversDecember 8, 2020 35/39 . Time step size series, one dimensional (cid:107) error (cid:107) (cm − ) convergence rate∆ t (s) L L ∞ L L ∞ . . − . − . . . .
210 0 . . . .
16 1 . . . .
28 2 . .
120 0 .
265 2 .
95 2 . .
930 2 .
02 3 .
08 3 . .
88 17 . B. Spatial point distance series, one dimensional (cid:107) error (cid:107) (cm − ) convergence rate∆ x (cm) L L ∞ L L ∞ . . .
216 0 . . . .
28 1 . . . .
43 3 . .
263 0 .
516 3 .
95 3 . .
07 7 . Table S1.
Convergence test results, one spatial dimensionalEqs. (53, 11, 12) were solved numerically from t = 0 s to 8192 s on x ∈ Ω = [0 , δρ (45) grew from a ρ = 1 to a ρ e λ ≈ A shows the results of varying the time step size from 4 to 256 s (with a fixed spatial pointdistance of ∆ x =1 /
512 cm). B shows the results of varying the spatial point distancefrom 1 / /
64 cm (with a fixed time step of 4 s). The error in ρ at the final timepoint was calculated as the difference between the numerical result and exact result (46). L and L ∞ norms of the error are tabulated. Convergence rate is calculated betweenconsecutive rows as log( (cid:107) error (cid:107) / (cid:107) error (cid:107) ) / log( h /h ), with h being either ∆ t or ∆ x ,as appropriate. The mean of ρ was 9000 cm − in all cases. Thus the relative erroris about 1 / . / ≈ . × − for ∆ t = 4 s,∆ x = 1 /
512 cm in one dimension. Errors in U a and U r (not shown) were smaller butotherwise behaved similarly. All numerical solutions used the PETSc Rosenbrock-Wsolver ra34pw2 (nominally an order 3 method), and fourth-order approximations for thespatial derivatives.December 8, 2020 36/39 . Time step size series, two dimensional (cid:107) error (cid:107) (cm − ) convergence rate∆ t (s) L L ∞ L L ∞ . . .
01 2 .
208 0 . .
215 0 .
692 0 . . . − . − . . . .
40 3 . .
167 0 .
800 1 .
53 1 . .
483 1 .
94 3 .
14 3 . .
27 18 . B. Spatial point distance series, two dimensional (cid:107) error (cid:107) (cm − ) convergence rate∆ x (cm) L L ∞ L L ∞ √ / . . − . − . √ / . . . − . √ /
512 0 . . .
84 2 . √ /
256 0 . .
199 3 .
92 3 . √ /
128 0 .
838 2 . Table S2.
Convergence test results, two spatial dimensionsEqs. (53, 11, 12) were solved numerically from t = 0 s to 8192 s on ( x, y ) ∈ Ω =[0 , √ / × [0 , / δρ (49) grew from a ρ = 1to a ρ e λ ≈ A shows the results of varying the time step size from 4 to 256s (with a fixed spatial point distance of ∆ x = √ / B shows the results ofvarying the spatial point distance from √ / √ /
128 cm (with a fixed time stepof 4 s). In all cases ∆ y = ∆ x × √ / igure S1. Attractant-only simulation reruns A. These ten images reproduce the numerical experiment of Figure 2B—simulation ofthe attractant-only model in two dimensions—but with different pseudorandom noisein the initial condition. Only worm density ρ ( x, y ) at t = 200 000 s (2 days, 7:33:20) isshown. B. Like A , but at t = 1 × s (115 days, 17:46:40). These images correspondone-to-one to the images in A . The number of aggregates in these panels ranges from oneto three, although a single aggregate may appear in as many as four pieces because ofthe periodic boundary conditions. To ease the identification of aggregates, the aggregateto which each piece belongs is identified by a red number.720 min. There is one frame per minute of real time, and the playback rate is7 s − . options138a.mp4 Numerical solution of attractant-only model in one dimension.This video corresponds to Fig 2A. This and all following videos are 200s long at15 s − . Time is displayed as “days, H:MM::SS”. Time ranges from 0 s to 200 000 s(2 days, 7:33:20). The two numbers below each panel are the minimum andmaximum of the plotted field. options139.mp4 Numerical solution of attractant-only model in two dimensions.December 8, 2020 38/39 igure S2.
Attractant+repellent simulation reruns A. These ten images reproduce the numerical experiment of Fig 3B—simulation of theattractant+repellent model in two dimensions—but with different pseudorandom noisein the initial condition. Only worm density ρ ( x, y ) at t = 200 000 s (2 days, 7:33:20) isshown. B. Like A , but at t = 1 × s (115 days, 17:46:40). These images correspondone-to-one to the images in A .Corresponds to Fig 2B. options140a.mp4 Numerical solution of attractant+repellent model in one dimension.Corresponds to Fig 3A. options141.mp4
Numerical solution of attractant+repellent model in two dimensions.Corresponds to Fig 3B. options157.mp4
Numerical solution of attractant+repellent model in two dimensionson a 6 cm ××