Analysis and Simulation of a Novel Run-and-Tumble Model with Autochemotaxis
BBulletin of Mathematical Biology Preprint
A Run-and-Tumble Model with Autochemotaxis
Nicholas J. Russell and Louis F. Rossi
September 5, 2020
Abstract
A specific species of phytoplankton,
Heterosigma akashiwo , hasbeen the cause of harmful algal blooms (HABs) in waterways around theworld causing millions of dollars in damage to farmed animals and destroyingecosystems. Developing a fundamental understanding of their movements andinteractions through phototaxis and chemotaxis is vital to comprehending whythese HABs start to form and how they can be prevented. We develop a one-and two-dimensional mathematical and computational model reflecting themovement of an ecology of plankton, incorporating both run-and-tumble mo-tion and autochemotaxis. We present a succession of complex and biologicallymeaningful models combined with a sequence of laboratory and computationalexperiments that inform the ideas underlying the model. By analyzing the dy-namics and pattern formation which are similar to experimental observations,we identify parameters that are significant in plankton’s pattern formation inthe absence of bulk fluid flow. We numerically analyze variations on how mightplankton deposit chemical and connect the outcomes with features observedin experimental observations.
Keywords
Plankton · Run-and-tumble · Chemotaxis · Heterosigma akashiwo
Plankton is a general term for a wide range of floating or drifting organ-isms in fresh and salt water. They are not strong enough to swim againstcurrents. Plankton are of considerable scientific interest for a number of rea-sons. For one, they are an important food source for larger organisms. Also,
Nicholas J. Russell · Louis F. RossiDepartment of Mathematical SciencesUniversity of DelawareNewark, DE 19716Tel.: (302) 831-4446E-mail: [email protected] a r X i v : . [ m a t h . D S ] S e p Nicholas J. Russell and Louis F. Rossi they process about 60% of the marine carbon, and they consume a significantfraction of man-made CO2 while producing about half of the atmosphericoxygen (Berdalet et al. 2016; Riebesell 2008). Nonetheless, plankton exhibita rich array of behaviors including self-propulsion, aggregation, photosynthe-sis, feeding, and predation (Sheng et al. 2007, 2010). Sometimes, coordinatedmicroscopic behaviors lead to significant mesoscale processes. Most notably,several types of plankton are responsible for harmful algal blooms (HABs)wherein the population density of plankton skyrockets in a concentrated area(Hallegraeff et al. 2004). However, some types of plankton, and particularlythose involved in the formation of HABs, aggregate even when there is noexternal current. While it is reasonable to model plankton as independentpassive scalar quantities drifting in a flow field, a more detailed understandingof their behavior may help scientists understand and anticipate when planktonwill aggregate and what combinations of features will drive this process.In this paper, we model, analyze, and simulate an autochemotactic plank-ton popoulation undergoing run-and-tumble swimming using the phytoplank-ton
Heterosigma akashiwo as our model organism. Over the past several years,this raphidophyte has caused several HABs throughout the world, includingJapan, Scandanavia, South Carolina, and the Delaware Bay (Hennige et al.2013). These HABs are responsible for large fish kills (e.g. the Puget Sound in2006) and public health hazards for the communities that surround the coasts(Lewitus et al. 2012). This is of particular interest to the Delaware commu-nity, as these HABs are becoming more frequent over the past several yearsin the Delaware Bay.
Heterosigma akashiwo are predatory swimmers, capableof consuming bacteria and performing an inefficient photosynthesis. Even inthe absence of prey or bulk fluid motion, plankton exhibit complex patternformations in two- and three-dimensions. Their behavior is a direct responseof local chemical and optical cues, so a fundamental understanding of howthey swim and hunt is essential to anticipating changes in marine ecology inresponse to climate change and ocean acidification (Fu et al. 2012).In this paper, we discuss a research program at the University of Delawareto mathematically model, analyze, and understand complex plankton behav-ior. These activities include careful experiments, isolating specific variables andbehaviors to connect them to a unique coupling of well-established processesin predation, chemotaxis, and phototaxis. We develop and analyze a model forplankton motion and signaling that is built upon the works of several authors(Hillen 2002; Hillen and Painter 2009; Keller and Segel 1971a,b; Lushi et al.2012; Bearon and Pedley 2000) in order to gain insight and understanding intothe interacting mechanisms of organism propulsion and chemical deposition,diffusion, and decay.The motion of plankton is driven by a number of processes, only some ofwhich have been fully explored experimentally for our model organism. There-fore, it is necessary to limit the study of their motion to a subset of these whichhave the greatest impact on their bulk coordinated behavior. In this paper,we will limit our discussion to run-and-tumble propulsive behavior which hasbeen studied extensively (Keller and Segel 1971a,b). When combined with au-
Run-and-Tumble Model with Autochemotaxis 3 tochemotaxis, in which organisms release the chemical they are attracted to,run-and-tumble motion admits quite rich dynamics which mirror experimentalfindings. We develop and analyze a one- and two-dimensional model of plank-ton motion incorporating run-and-tumble motion along with autochemotaxiswhile using a biologically relevant tumbling probability.In Section 2, we summarize the experiments conducted and the resultantobservations which inform our mathematical model. In Section 3, we developanalytical models of run-and-tumble dynamics coupled with autochemotaxisin one and two dimensions. In Section 4, an analysis of the one-dimensionalmodel is conducted both analytically and numerically. In Section 5, we ana-lyze our two-dimensional model analytically and numerically. In Section 6, wediscuss the analysis and simulations by connecting our results to experimentalobservations. Finally, in Section 7, we share our conclusions of this work andfuture research opportunities.
In this paper, we use
Heterosigma akashiwo , a constituent plankton in an algalbloom colloquially referred to as the “red tide”, as our model organism. Togain insight on the evolution and dynamics of the plankton distributions, weperformed a series of laboratory experiments using a camera to gather rawdata. See Appendix A for the full details on experimentation. The data shouldbe considered qualitative in the sense that we could not rigorously map theimage density to a calibrated population density.We anticipate that two main mechanisms drive the plankton density dis-tribution in water without external fluid flow: chemotaxis and phototaxis. Forchemotaxis, plankton swim in response to a chemical while also emitting thechemical, also called autochemotaxis (Seymour et al. 2009). Phototaxis is de-pendent on the plankton’s photosystem behavior and the depth of the water,and light is a source of energy for the organism. Its photosystem is degradedby light and is constantly repaired by internal biochemical reactions (Fu et al.2008). As well, plankton are observed to vertically and horizontally migratetowards and away from light, and it is reasonable to hypothesize that thestate of their photosystem drives this behavior (Hara and Chihara 1987). Itis our aim to investigate autochemotactic behavior in this model organism asa foundation for exploring its richer behavior in response to phototaxis andother processes impacting migration and aggregation.To demonstrate the richness of the full dynamics of
Heterosigma akashiwo ,even in the absence of an external driving fluid flow, we conducted experimentsin a deep tank. In this deep tank, all effects from chemotaxis and phototaxisare observable. We observe plankton aggregating into complex patches andfilaments, while migrating towards and away from the light (see Fig. 1a).Water attenuates light, captured by the Lambert-Beer Law, and aggregationscan shade other plankton.
Nicholas J. Russell and Louis F. Rossi (a)
Deep plate experiment (b)
Shallow plate experiment
Fig. 1
Pictures from various experiments conducted. The aggregations in b formed after 5minutes in complete darkness (an infrared camera was used). To remove optical attenuation and shading from consideration and to min-imize the impact of vertical migration, the depth of the water column is re-duced until the vertical effects were no longer observable, which occurs whenthe depth is less than 1 cm. All experiments we discuss from now on onlyconsider a shallow dish to remove vertical migration from consideration as adominant affect. The system is flooded with light, and thus there is no photo-intensity gradient in the field. Therefore, the shallow environment allows us toeliminate certain affects attributable to the organisms sensitivity to light. Westill must account for the fact that
Heterosigma akashiwo responds to light,even if light is uniformly distributed in the shallow environment.To further isolate the behavioral impact of photosynthesis, we conduct anexperiment where we film
Heterosigma akashiwo in complete darkness, allow-ing the photosystems in all organisms to reset, and then observe what occursonce light returns. Within the first 10 minutes, aggregations form in completedarkness (see Fig. 1b). Thus, we can establish that autochemotaxis drives theformation of aggregations. Conducting experiments with lights on, we observethat phototaxis drives the coarsening and movement of these aggregations.From a simple experiment shown in Fig. 2, aggregations started forming within2 minutes, more rapidly than the time it takes the photosystems to degradeand break. The structure of the aggregations is also different with and withoutlight, so phototaxis changes the dynamics of aggregation shapes.We conducted experiments with a lid flush with the upper liquid surfaceto see if Marangoni effects were playing a role in the formation of aggrega-tions and found no observed difference with all other free surface experiments.From these experiments and observations, we infer that although phototaxis
Run-and-Tumble Model with Autochemotaxis 5 (a)
Time Elapsed: 2 Minutes (b)
Time Elapsed: 6 Minutes (c)
Time Elapsed: 15 Minutes (d)
Time Elapsed: 40 Minutes (e)
Time Elapsed: 60 Minutes
Fig. 2
The evolution of aggregationpatterns from
Heterosigma akashiwo during an hour long experiment. Weused a square, shallow dish with alength of 6.3 cm along with a lampplaced at a 60 ◦ angle of incidence.We used a filter on each picture tomake the aggregations clearer. Wenote that the spiral pattern in a iscaused by initial fluid motion whenthe plankton solution is poured intothe dish initially. Nicholas J. Russell and Louis F. Rossi stimulates motion and perhaps aggregation, chemotaxis combined with run-and-tumble is the initial driver of aggregations of Heterosigma akashiwo . N p plankton as point sources, ρ ( x , t ) = N p (cid:88) i =1 δ ( x − x i ( t )) , (1)where x i ( t ) = ( x ( t ) , x ( t )) is the position of the i th plankton and δ ( x ) isthe Dirac delta function. We denote c ( x , t ) as the chemical concentration ata point x . When considering a classic run-and-tumble model for planktonmotion, literature recommends the flipping probability of P [Tumble in ( t, t + ∆ t )] = λ (cid:18) − v · ∇ c (cid:107) v · ∇ c (cid:107) (cid:19) ∆ t (2)where v is the velocity of a plankton, λ ≥
0, and ∇ c is the chemical gradient(Saintillan 2018). This model is problematic from a biological perspective. Wesee that v ·∇ c (cid:107) v ·∇ c (cid:107) = sgn ( (cid:107) v · ∇ c (cid:107) ) takes three distinct values: − , , and 1. Thisgives a drastic shift in behavior when plankton flips to being orthogonal tothe gradient, implying an infinite sensitivity to both the environment and theorganism’s response to the chemical. Thus, we propose to create a transitionto smooth out this hard shift by using the parameter δ > P [Tumble in ( t, t + ∆ t )] = λ − v · ∇ c (cid:113) ( v · ∇ c ) + δ ∆ t. (3)Biologically, δ captures the ability of an organism to sense weak gradients inits environment. As δ →
0, the organism can perfectly sense an infinitesimallysmall gradient, and we recover our original model (2). Larger values of δ de-scribe a smooth transition in the organism’s behavior which could capture anoisier response to weak signals when averaged over multiple instances of thesame organism or multiple organisms in close proximity. We shall see laterthat δ has a strong impact on macroscopic pattern formation.We assume the plankton move in a run-and-tumble motion with chemo-taxis, i.e. the probability of a plankton tumbling and changing directions isgiven by (3). If a plankton does tumble, the probability that a plankton movesfrom a direction v to v (cid:48) is defined by the transition probability T ( v (cid:48) , v ). Weassume that plankton move at a constant speed v and thus, v = v e θ , where e θ = (cid:104) cos( θ ) , sin( θ ) (cid:105) . Because of this assumption, the transition probabilitycan be rewritten as T ( v (cid:48) , v ) = T ( θ (cid:48) , θ ). Run-and-Tumble Model with Autochemotaxis 7
The chemical attractant, which is produced by each plankton, diffuses anddecays throughout the system. Since there is no consensus on when the plank-ton emit this chemical, we will assume that plankton emit this chemical basedon how much chemical they can sense at their current position. Let f ( c ) be ageneral deposition function that is only dependent on the chemical, c . The dif-fusion and decay of c can be captured in an evolution equation for the chemicalattractant: (4) c t = κ ∆ c − βc + f ( c ) ρ, where κ ≥ β ≥ ρ from an Eulerian perspective, we now define ρ ( x , t ) to be the density of the plankton at a given time and location in R .We introduce ψ ( x , θ, t ) as the density of plankton at time t and position x ,moving in the θ direction. Thus, the evolution equation for ψ can be writtenas a Fokker-Plank equation: (5) ρ ( x , t ) = (cid:90) π ψ ( x , θ, t ) dθ (6) ψ t = − v e θ · ∇ ψ − λ − v ( ∇ c · e θ ) (cid:113) v ( ∇ c · e θ ) + δ ψ + λ (cid:90) π − v ( ∇ c · e θ (cid:48) ) (cid:113) v ( ∇ c · e θ (cid:48) ) + δ T ( θ, θ (cid:48) ) ψ ( x , θ (cid:48) , t ) dθ (cid:48) . In (6), the first term represents plankton that have moved away from position x ; the second term represents the plankton at x who tumbled to a differ-ent direction than θ ; and the third term represents the plankton that havetumbled at x to the direction of θ from θ (cid:48) . Therefore, (4)-(6) is the contin-uum two-dimensional run-and-tumble with autochemotaxis system. To non-dimensionalize this model, we denote the natural time-scale of τ = λ and thenatural length scale of L = vλ . Define the non-dimensional time and length as t = τ t (cid:63) , x = L x (cid:63) , where the (cid:63) -notation denotes a non-dimensional parameter.To simplify further, we assume that the turning kernel is a uniform distribu-tion, i.e. T ( θ (cid:48) , θ ) = π . Removing the star notation, we obtain the followingnon-dimensional model of autochemotactic-driven motion in 2D: (7a) c t = d ∆ c − d c + f ( c ) ρ (7b) ρ ( x , t ) = (cid:90) π ψ ( x , θ, t ) dθ (7c) ψ t = − e θ · ∇ ψ − − ∇ c · e θ (cid:113) ( ∇ c · e θ ) + δ ψ + 14 π (cid:90) π − ∇ c · e θ (cid:48) (cid:113) ( ∇ c · e θ (cid:48) ) + δ ψ ( x , θ (cid:48) , t ) dθ (cid:48) , Nicholas J. Russell and Louis F. Rossi where d = κλ v , d = βλ , and f ( c ) has been scaled by λ . Throughout thispaper, we will analyze and simulate the system in a large, square domain[0 , (cid:96) ] × [0 , (cid:96) ] with periodic boundary conditions to express the near limitlessdomain in which plankton truly live. Computational limitations necessarilyrequire that we make (cid:96) finite.3.2 The 1D modelFor a one-dimensional model, θ is limited to values of 0 (right-moving) and π (left-moving), and ψ , ρ , and c vary with x ∈ Ω ⊂ R and t ≥ θ dependence entirely by replacing ψ ( θ, x, t ) with the distinctfunctions ψ + ( x, t ) and ψ − ( x, t ), representing right- and left-moving planktonrespectively. Equation (3) can be rewritten as P [Tumble in ( t, t + ∆ t )] = λ (cid:32) − c x (cid:112) c x + δ (cid:33) , (8)and we assume that the plankton move at a constant speed v . For the au-tochemotaxis part of the model, let f ( c ) denote an arbitrary deposition func-tion. We can then write the evolution equations as follows: (9a) c t = κc xx − βc + f ( c )( ψ + + ψ − ) (9b) ψ + t = − vψ + x − λ (cid:32) − c x (cid:112) c x + δ (cid:33) ψ + + λ (cid:32) c x (cid:112) c x + δ (cid:33) ψ − (9c) ψ − t = vψ − x + λ (cid:32) − c x (cid:112) c x + δ (cid:33) ψ + − λ (cid:32) c x (cid:112) c x + δ (cid:33) ψ − . Non-dimensionalizing, we denote the natural time-scale of τ = λ and thenatural length scale of L = vλ . Define the non-dimensional time and lengthas t = τ t (cid:63) and x = Lx (cid:63) , where the (cid:63) -notation denotes a non-dimensional pa-rameter. Removing the star notation, we obtain the following non-dimensionalmodel of (9a) − (9c): (10a) c t = d c xx − d c + f ( c )( ψ + + ψ − ) (10b) ψ + t = − ψ + x − (cid:32) − c x (cid:112) c x + δ (cid:33) ψ + + 12 (cid:32) c x (cid:112) c x + δ (cid:33) ψ − (10c) ψ − t = ψ − x + 12 (cid:32) − c x (cid:112) c x + δ (cid:33) ψ + − (cid:32) c x (cid:112) c x + δ (cid:33) ψ − where d = κλ v and d = βλ , and f ( c ), which is scaled by λ , are all non-dimensional. If we define ρ ( x, t ) = ψ + + ψ − to be the total density, then Run-and-Tumble Model with Autochemotaxis 9
Fig. 3
Examples of the three types of deposition functions considered in the simulations.The black dashed line is f ( c ) = γ , the constant deposition function; the green, solid line is f ( c ), the switch deposition function; and the orange, dotted line is f ( c ), the linear switchdeposition functions. The analytic descriptions of f ( c ) and f ( c ) can be found in equations(12a) and (12b), respectively. manipulating (10b) and (10c) allows us to obtain a coupled PDE system for ρ and c : (11a) c t = d c xx − d c + f ( c ) ρ (11b) ρ tt + ρ t = ρ xx − ∂∂x (cid:34) c x (cid:112) c x + δ ρ (cid:35) For both (10a)-(10c) and (11a)-(11b), we will assume periodic boundary con-ditions and consider the domain Ω = [0 , (cid:96) ], where (cid:96) >
Heterosigma akashiwo is unknown, wewill utilize basic deposition functions related to other organisms such as insectsand bacteria. The three deposition functions we analyze are f ( c ) = γ , f ( c ) = γ Θ( c − c ), and f ( c ) = γ ( c + η )Θ( c − c ), where γ, c , η > x ) is theHeaviside function. The constant function, f , is where an organism alwaysemits the chemical; the switch, f , is where an organism releases the chemicalonly if the chemical concentration at x is below the threshold c ; and thelinear switch, f , is where an organism releases a specific amount of chemicaldependent on the amount of chemical at x , but does not emit if the chemicalat x is above c . In our numerical simulations, we smooth out the Heavisidefunctions to avoid discontinuities as (12a) f ( c ) = γ (cid:20) tanh (cid:18) c − cT (cid:19) + 1 (cid:21) (12b) f ( c ) = γ ( c + 0 . c )2 c (cid:20) tanh (cid:18) c − cT (cid:19) + 1 (cid:21) , κ Chem. diffusion rate 10 − − − cm / s Stocker and Seymour 2012 β Chem. decay rate 5 −
40 s − Kawaguchi 2011 λ Tumbling rate 2 −
50 s − Visser and Thygesen 2003 v Swimming speed 1 × − cm/s Harvey et al. 2015 Table 1
Dimensional parameters with descriptions and experimental values, along withcitations.Parameter Description Important Details Simulation Values d Chem. diffusion rate κλ v . − d Chem. decay rate βλ . − δ Run-and-Tumble parameter – 0 − . (cid:96) Length of domain – 3 − γ Max. chemical strength – 0 . − . c Chem. sensitivity threshold f and f only 0 . − . T Behavior switch length – 0 . − . N p Number of plankton 2D system only 10 − Table 2
Non-dimensional parameters used in 1D and 2D simulations along with descrip-tions, expressions of values in terms of dimensional parameters, and important details, alongwith the numerical values used in our simulations. where T >
0. Examples of these deposition function are displayed in Fig. 3.For a summary of all dimensional and non-dimensional parameters, see Tables1 and 2, respectively. − (10c). There is a constant steady state from(10a) of d c = f ( c ) ρ , where ρ = 2 ψ and ψ is the constant solution for both ψ + and ψ − . We linearize around this solution by letting c = c + (cid:15)b, ψ + = ψ + (cid:15)a + , ψ − = ψ + (cid:15)a − , (13)where (cid:15) (cid:28) a + , a − , and b are functions in terms of time and space.Substituting (13) into (10a) − (10c), at leading order of (cid:15) we obtain (14a) b t = d b xx + d b + f ( c )( a + + a − ) (14b) a + t = − a + x + 12 ( a − − a + ) + ψb x δ (14c) a − t = a − x −
12 ( a − − a + ) − ψb x δ , Run-and-Tumble Model with Autochemotaxis 11 where d = ρf (cid:48) ( c ) − d . As a note, we utilized the Taylor expansion x √ x + δ = xδ − x δ + O ( x ). We now consider the Fourier transform of the (14a)-(14c),and let A + , A − , and B denote the transforms of a + , a − and b , respectively.Letting k denote the wave number in the Fourier transform, we obtain (15a) B (cid:48) = − k d B + d B + f ( c )( A + + A − ) (15b)( A + ) (cid:48) = − ikA + + 12 ( A − − A + ) + ikψδ B (15c)( A − ) (cid:48) = ikA + −
12 ( A − − A + ) − ikψδ B, which can be rewritten into the matrix form BA + A − (cid:48) = d − d k f ( c ) f ( c ) ikψδ − ik −
12 12 − ikψδ ik − BA + A − . (16)We now look for eigenvalues of the 3 × λ as aneigenvalue of the above matrix, we can write the characteristic equation as(17) λ + (cid:2) d k − d + 1 (cid:3) λ + (cid:2) ( d + 1) k − d (cid:3) λ + k (cid:18) d k − d − cd δ (cid:19) = 0 . We seek to find combinations of parameters in which the constant solutionis unstable. Take the roots of (17) as the set { λ ( k ) , λ ( k ) , λ ( k ) } and define R ( k ) := max { Re ( λ i ) } i =1 . Let k u be the most unstable wave number, i.e. k u = arg max k> R ( k ) . (18)If k u > π(cid:96) , the set of parameters will make ( c, ρ ) unstable in the domain length (cid:96) . Given a fixed set of parameters d , d , c, δ and deposition function f ( c ), wecan vary k to numerically find the most unstable wave number. Examples ofthe function R ( k ) are shown in Fig. 4. In Fig. 4a, we keep d constant but varythe diffusion rate, d . As d increases, k u decreases. Alternatively, in Fig. 4b, d is kept constant but the chemical decay rate, d , is varied. As d increases, k u increases.To see this relationship between d and d more clearly, we plot the stabilityregions for a domain length (cid:96) in the d - d plane, shown in Figs. 4c and 4d.The green, shaded regions are examples of parameter regimes where the systemwould be unstable when (cid:96) = 6 for either δ = 0 .
01 or δ = 0 . δ increasesincrementally, the stable regime becomes much larger. As well, we note thatthe curve separating the unstable and stable regimes becomes linear as d and d both increase.We now remark on some properties of (17) and R ( k ), the proofs of whichcan be found in Appendices B.1 and B.2, respectively. These two propositionsgive stability conditions for large wave numbers and when the constant solution( c, ρ ) may admit an unstable solution. (a) Varying d with d = 1 (b) Varying d with d = 1 (c) Stability regions with δ = 0 . (d) Stability regions with δ = 0 . Fig. 4
Plots a and b display R ( k ) with varied d and d parameters. Plots c and d showphase space stability for a constant solution of (11a)-(11b) with varied δ . Parameters usedcan be found in Appendix C.2. Proposition 1 If d > , lim k →∞ R ( k ) = − . Proposition 2
Define d = ρf (cid:48) ( c ) − d . Then, if − cd δ < d < d (cid:34) (cid:114) cd d δ (cid:35) , (19) there will be an (cid:96) such that the constant solution ( c, ρ ) to (10a) − (10c) isunstable. Run-and-Tumble Model with Autochemotaxis 13
Fig. 5
Time-series of plots for simulation of 1D autochemotactic system (11a)-(11b) at t = 0 , , , ρ . The middle rowdisplay the evolution of chemical concentration, c . The bottom row shows the evolution of E ( k ) (see (20)), where k are the Fourier modes. The black vertical lines in the bottom roware at the most unstable wave number, k u ≈
8. This simulation can be seen in its entirelyin Online Resource 1 and parameters used can be found in Appendix C.2.
In Figure 5 and Online Resource 1, we observe the evolution of the systemin a parameter regime where the constant stationary solution is unstable. Theunstable distribution of aggregations rapidly saturates from nonlinear inter-actions and yields a slowly evolving coarsening pattern, eventually becominga single aggregation. The top two panels display the plankton density, ρ , andchemical density, c , respectively. The lower panel displays E ( k ), which showsthe real part of the Fourier modes of c , i.e. E ( k ; t ) = Re {F [ c ( x, t )] ( k ) } , (20)where F is the Fourier transform operator. To calculate this numerically, weutilize a fast Fourier transform. Aside from the translational mode of k = 0,the most unstable wave number, k u , will grow the quickest initially. In thissimulation, the most unstable wave number is k u ≈
8, denoted by the black,dashed lines on the plots. At t = 4, six small aggregations are formed and themode at k = k u has been triggered, as it is the most pronounced k at the outset.At t = 20, the aggregations have coarsened down to only 2 aggregations, and k u has diminished from view due to this clustering. The final panel at t = 150shows that the two aggregations have merged, forming one large aggregation.We emphasize the timescale difference between the coarsening from six to twoaggregations and from two to one aggregations. (a) Simulation using f ( c ) = f ( c ) (animation in Online Resource 1) (b) Simulation using f ( c ) = f ( c ) (animation in Online Resource 2) (c) Simulation using f ( c ) = f ( c ) (animation in Online Resource 3) (d) Total chemical over time, C ( t ) Fig. 6
Time-series of plots for simulation of 1D autochemotactic system (11a)-(11b) forvarious deposition functions. All simulation parameters, which can be found in AppendixC.2, result in an unstable constant solution with k u = 8 . a . The total chemical over time for all three simulationsis in d . In Figs. 6 and 7, we explore two important facets of our model: the depo-sition function and δ . Recall our three deposition functions defined in (12a)and (12b). For this comparison, we choose parameter values for f ( c ), f ( c ),and f ( c ) such that all systems are unstable and have the same k u .Some key differences emerge from these simulations in Figs. 6a-6c. Eventhough all simulations have a similar number of aggregations from t = 4 to Run-and-Tumble Model with Autochemotaxis 15 t = 20, ending up with two aggregations at t = 20, the time needed to fullymerge the final two aggregations varies. The linear switch deposition function, f , admits a singular aggregation at t = 100, but it takes until t = 350 for theswitch deposition function, f , to merge. As well, we see a clear difference in theprofile for c achieved at the non-constant steady state. The f function admitsa profile for c and ρ that are relatively similar, but c and ρ look vastly differentfor the f and f cases. This can be attributed to the threshold parameter, c , hampering the plankton’s ability to deposit chemical once c becomes toolarge.Other insights can be gleaned from the total chemical concentration profile C ( t ) (Fig. 6d). The total chemical for the f simulation remains at its equi-librium value. Between merging events for the f and f simulations, there issteady chemical total; however, several time steps prior to a merging event,there is a drastic decrease in the amount of chemical due to the restructur-ing of ρ . We show that C ( t ) stays constant for f ( c ) = f ( c ) in a generaltwo-dimensional case in Proposition 3 in Section 5.This model is also highly sensitive to δ , which expresses the ability of theorganism to admit strong differentiated responses to small gradients. In Fig.7, the δ parameter is decreased to show the differences in the final steady statesolution of the model. As δ decreases, the peak of the plankton density andchemical concentration become more sharp. At δ = 0, note that a steady statesolution of (11a)-(11b) must satisfy ρ xx − [sgn( c x ) ρ ] x = 0 , (21)with periodic boundary conditions. Due to the translational invariance of thesolution to (21), we can center the aggregation in the middle of our domain[0 , (cid:96) ]. Using Fig. 7a, we recognize that sgn( c x ) > ≤ x ≤ (cid:96)/ c x ) < (cid:96)/ ≤ x ≤ (cid:96) . Thus, we obtain a steady state solution ofstitched exponential functions: ρ ( x ) = (cid:40) A e x , x ∈ [0 , (cid:96)/ A e − ( x − (cid:96) ) , x ∈ [ (cid:96)/ , (cid:96) ] , A = ρ(cid:96) (cid:0) e (cid:96)/ − (cid:1) . (22)The stability of this steady state is still unknown, along with the other non-constant steady state solutions with δ >
0. However, through numerical ex-perimentation, we conjecture that this steady state is unconditionally stable. d c = f ( c ) ρ and ρ = 2 πψ , (a) Plankton density with varied δ (b) Chemical concentration with varied δ Fig. 7
The non-constant steady states for the chemical concentration and plankton densityfor varied δ = 0 , . , . , . similar to the one-dimensional analysis. We perturb this constant solution as ψ = ψ + (cid:15)φ, c = c + (cid:15)µ, ρ = ρ + (cid:15) (cid:90) π φ := ρ + (cid:15)τ, (23)where φ, µ, τ are functions of time and space and (cid:15) (cid:28)
1. Substituting theperturbations from (23) to (7a), we obtain the evolution equation for µ at firstorder as µ t = d ∆ µ + d µ + f ( c ) τ, (24)where d = f (cid:48) ( c ) ρ − d . By substituting the perturbations into (7c) and uti-lizing the same Taylor expansion as in the one-dimensional case for the termswith δ , we obtain (25) (cid:15)φ t = − (cid:15) e θ · ∇ φ − (cid:18) − (cid:15) ∇ µ · e θ δ (cid:19) ( ψ + (cid:15)φ )+ 14 π (cid:90) π (cid:20) − (cid:15) ∇ µ · e θ (cid:48) δ (cid:21) ( ψ + (cid:15)φ ( θ (cid:48) )) dθ (cid:48) . At first order, we simplify (25) to become φ t = − e θ · ∇ φ − φ + ψ δ ∇ µ · e θ + 14 π τ. (26)Taking the two-dimensional Fourier transform in both x and y of (24) and(26) while denoting k = (cid:104) k , k (cid:105) as the wave vector, we get (27a) H t = i ( k · e θ ) H − H − ψi δ ( k · e θ ) M + 14 π T (27b) M t = ( d − d | k | ) M + f ( c ) T, Run-and-Tumble Model with Autochemotaxis 17 where H , M , and T are the two-dimensional Fourier transforms of φ , µ , and τ , respectively. Continuing, we write H in terms of its Fourier coefficients withrespect to θ , i.e. H ( k , t, θ ) = (cid:88) n ∈ Z h n ( k , t ) e inθ , (28)where h n ∈ C . From this representation, T can be written as T ( k , t ) = (cid:90) π H ( k , t, θ (cid:48) ) dθ (cid:48) = 2 πh ( k , t ) . (29)We define the complexification of the wave number k as ω := k + ik , anddenote ω as the complex conjugate of ω . We rewrite sin( θ ) and cos( θ ) inexponential forms and substitute (28) and (29) into (27a)-(27b). By matchingpowers of e inθ , we obtain the following infinite linear ODE system, whereprimes denote derivative with respect to t : (30a) h (cid:48) = i ωh − + ωh ) (30b) h (cid:48) = i ωh + ωh ) − h − ψiω δ M (30c) h (cid:48)− = i ωh − + ωh ) − h − − ψiω δ M (30d) h (cid:48) j = i ωh j − + ωh j +1 ) − h j (30e) M (cid:48) = ( d − d | k | ) M + 2 πf ( c ) h , where j (cid:54) = − , ,
1. This system can be rewritten as matrix equation withinfinite entries by y (cid:48) = Ay , where y = ... h − h − h Mh h ... , A = . . . . . . . . . i ω − i ω i ω − i ω − ψiω δi ω i ω πf ( c ) d − d | k | i ω − ψiω δ − i ω i ω − i ω . . . . . . . . . (31)We now look for the eigenvalues of A which have the largest real part. Sincewe cannot compute every eigenvalue of A algebraically, we will compute themnumerically by truncating the matrix at an appropriate threshold. Denote A N and y N as the truncated matrix A and vector y in (31) consisting of theequations with the functions { h − N , · · · h − , h , M, h , · · · , h N } , where N ∈ Z + is sufficiently large (see next paragraph). We now calculate the eigenvalues ofthe system y (cid:48) N = A N y N . Denote the set of these eigenvalues as { λ i ( | k | ) } N +2 i =1 , and define R N = R N ( | k | ) := max { Re ( λ i ) } N +2 i =1 . Then, K u will be denoted asthe norm of the most unstable wave number, i.e. K u := (cid:12)(cid:12)(cid:12)(cid:12) arg max k R N ( | k | ) (cid:12)(cid:12)(cid:12)(cid:12) . (32)If K u ≥ π(cid:96) , the set of parameters will make the constant solution ( c, ρ ) un-stable in the domain [0 , (cid:96) ] × [0 , (cid:96) ]. Given a fixed set of parameters d , d , c, δ and deposition function f ( c ), we vary k to find K u .Examples of the function R N ( | k | ) are shown in Figs. 8a and 8b, where theparameters d and d are varied. In a similar fashion to the one-dimensionalcase (Fig. 4), decreasing d or increasing d causes K u to increase and R N ( | k | )approaches − / | k | approaches infinity. In Fig. 8c, we found that N = 100fully resolved R N ( | k | ). This shows the rapid convergence of small | k | for all N and as N increases, the approximation of large | k | behavior converges towardsan asymptote at − /
2. We also highlight how c , the threshold parameterfor the deposition functions, changes the stability of the constant solution. InFig. 8d, by increasing c for f in (12a), the function R N ( | k | ) for a constantdeposition function with the same parameters is recovered.Finally, Figs. 8e and 8 f display the stability regions in the d - d plane fortwo different δ values: 0.001 and 0.004. As δ increases, there are fewer steadystates which will be unstable for a certain domain size [0 , (cid:96) ] × [0 , (cid:96) ]. As withthe one-dimensional case, the curve separating the stable and unstable regimesbecomes linear as d and d both increase.5.2 Linear Growth and Nonlinear SaturationWe now turn towards understanding the nonlinear saturation of unstablesteady state solutions to our two-dimensional model. We analyze the La-grangian system for ρ and c and show numerical results. Table 4 and Fig.3 show the parameters and three deposition function we will utilize in oursimulations, respectively.To initialize our system, the density ρ can be constructed as a summationof N p point sources, ρ ( x , t ) = N p (cid:88) i =1 δ ( x − x i ( t )) , (33)where x i is the position of each plankton and δ ( x ) is the Dirac delta function.Initially, all organisms are randomly distributed throughout the domain andgiven a random orientation θ i , selected from a uniformly random distribution[0 , π ). In each iteration, plankton tumble with a non-dimensional probabilityof P [Tumble in ( t, t + ∆ t )] = 12 − e θ i · ∇ c (cid:113) ( e θ i · ∇ c ) + δ ∆ t. (34) Run-and-Tumble Model with Autochemotaxis 19 (a)
Varying d with d = 1 (b) Varying d with d = 1 (c) Varying N with d = d = 1 (d) Varying c with d = d = 1 (e) Stability regions with δ = 0 . (f) Stability regions with δ = 0 . Fig. 8
Plots a - d display the function R N ( | k | ), varying d , d , N and c , respectively. For a , b , and d , we utilize N = 100. Plots e and f show phase space stability for a constantsolution of (7a)-(7c) with varied δ . Parameters used can be found in Appendix D.2. If a plankton tumbles, a random direction is then chosen utilizing the turningkernel T ( θ (cid:48) i , θ i ) = π . Each particle then updates their position by x i ( t + ∆ t ) = x i ( t ) + ∆ t (cid:104) cos( θ i ) , sin( θ i ) (cid:105) . (35)For autochemotaxis, we calculate f ( c ) at these new positions using ourdeposition function of choice. To deposit the chemical numerically, we place a Fig. 9
A time-series evolution of a 2D simulation described in Section 5.2 at times t =0 , . , , , and 20 for ρ , c , and (cid:101) E ( k ). Note that the scaling for (cid:101) E changes every panel. Thelightest shade is always 0, but the darkest shade is determined by the second largest valueof (cid:101) E ( k ) at each t . The black, dotted circle in the bottom row has a radius K u . For thissimulation, K u = 6 . two-dimensional Gaussian centered at x i instead of using a Dirac delta func-tion. We then evolve equation (7a) by utilizing a Crank-Nicolson time-steppingmethod. See Appendix D for more specific details regarding the simulations.In Fig. 9, we display the evolution of a parameter regime which makes theconstant solution ( c, ρ ) unstable and induce pattern formation. To visualize ρ more effectively, we portray all N p plankton as small, two-dimensional Gaus-sians centered at x i and plot the resulting density. The top row of Fig. 9 showsthis density evolution; the middle row displays the evolution of c , the chemicalconcentration; and the bottom row displays (cid:101) E ( k ), which calculates the realpart of the Fourier modes of c , i.e. (cid:101) E ( k ; t ) = Re {F [ c ( x , t )] ( k ) } , (36)where F is the two-dimensional Fourier transform. To calculate this numer-ically, we utilize a two-dimensional fast Fourier transform. Other than thetranslational mode of k = , a wave vector whose norm is equal to K u (see(32)) should grow the quickest initially. For this simulation, the norm of themost unstable wave number is K u ≈
6, denoted by the circle with radius K u in the lower panels.At t = 0 .
2, the development of small aggregations occur and several modesgrow quickly in the Fourier domain, the largest of which occurs at K u ≈ (cid:101) E ( k ) as wave vec-tors with smaller norms emerge, reflecting pattern formation. Unsurprisingly, Run-and-Tumble Model with Autochemotaxis 21 chemical concentrations are highest near the most densely populated areas ofthe plankton density.We now discuss the effects the various deposition functions have on thedevelopment and structure of aggregations. To equitably compare all threedeposition functions, parameters were found such that K u ≈ d , d , and δ the same (see Table 4). The evolution of the plankton densitiescan be seen in Figs. 10a-10c, where the color bar is centered in relation to ρ for all plots (Fig. 10d). At t = 3, there are more pronounced aggregations withthe constant deposition function, f , when compared with the other depositionfunctions. At t = 71, we see that both f and f , the linear switch, have morenoticeable aggregations rather than f , the switch. In the final frame when t = 150, the aggregations for f are the most defined whereas f admits theleast pronounced aggregations. For all simulations, merging events occur morefrequently early on rather than later, due to the aggregations’ sizes.We utilize two metrics to elucidate the differences between these depositionfunctions, the first of which is the total chemical in the system. Denote C ( t )as the total chemical, i.e. C ( t ) = (cid:90) (cid:90) D c ( x , t ) d x , (37)where D = [0 , (cid:96) ] × [0 , (cid:96) ] . In Fig. 10e, the evolution of C ( t ) is shown for the threesimulations. All simulations were initialized with c ( x ,
0) = c = 0 .
01 and thus, C (0) = c(cid:96) = 1. The long term behavior of all three simulations is decidedlyvaried, with f staying constant, f sharply declining before leveling out, and f steadily declining before leveling out. Proposition 3, the proof of which is inAppendix B.3, shows that the constant deposition function admits no changein chemical if the chemical is initialized at c . Proposition 3
Denote the total chemical concentration as C ( t ) , defined in (37) . Let d c = f ( c ) ρ be the constant steady state solution for (7a) - (7c) on D = [0 , (cid:96) ] × [0 , (cid:96) ] . Then, if f ( c ) = f ( c ) and C (0) = C where C ≥ , then lim t →∞ C ( t ) = c(cid:96) . (38)To quantitatively show differences between the structure of these depositionfunctions, we seek to understand the probability of finding a point in D witha given plankton density at a given time t , i.e. S t ( q ) = P [ ρ ( x , t ) = q | x i ∈ D ] . (39)If the system is initialized with ( c, ρ ) which satisfies d c = f ( c ) ρ , then S ( q ) = δ ( q − ρ ) , where δ is the Dirac delta function. In Fig. 10f, we graph S forall three deposition functions along with S for comparison. The three distri-butions show clear structural differences for long time behavior. To begin, f admits a skewed right distribution, showing that large extreme values of ρ arepresent and thus the aggregations are densely populated. In contrast, f and f are approximately normally distributed, but f has a much larger variation (a) Evolution of ρ with f ( c ) = f ( c ) (animation in Online Resource 4) (b) Evolution of ρ with f ( c ) = f ( c ) (animation in Online Resource 5) (c) Evolution of ρ with f ( c ) = f ( c ) (animation in Online Resource 6) (d) Colorbar for all simulations (e)
Total chemical concentration over time (f) S for varied deposition functions Fig. 10
Plots a , b , and c are time-series plots for ρ using several different depositionfunctions, using the colorbar shown in d , where ρ satisfies d c = f ( c ) ρ . Plot e shows thetotal chemical in the system over time. Plot f shows S for all deposition functions alongwith S ( q ) = δ ( q − ρ ) (see (39)). All simulation parameters, which can be found on Table 4in Appendix D.2, are constructed to have K u ≈
6. Run-and-Tumble Model with Autochemotaxis 23 (a)
Plankton density with varied δ at t = 150 (b) S for varied δ Fig. 11
Long time behavior ( t = 150) of varied δ = 10 − , − , − , − . a shows theplankton density and b shows S for all δ values. For parameters, see Table 4. than f . This implies that there are more densely populated aggregations forthe linear switch deposition function, whereas the switch deposition functiondistributes the plankton across aggregations more evenly.Finally, we see that δ has a profound impact on quasi-steady aggregationsafter the system relaxes. Plankton density plots at t = 150 are shown in Fig.11a. As δ →
0, the aggregations become more defined and there are morevisible peaks. We can see this more clearly by plotting S in Fig. 11b. As δ →
0, the probability density becomes more skewed right, implying that theaggregations are more densely populated when δ decreases. We now discuss our results from the previous sections by connecting themtowards the experimental observations in Section 2.Our stability analysis in Figs. 4 and Figs. 8 gives us insight of how a con-stant solution ( c, ρ ) can become unstable and initiate plankton aggregations.If we decrease d or increase d , we are able to make the system more unstable,i.e., as the chemical diffuses more slowly or decays more rapidly, the planktonare more likely to aggregate. Also, the importance of δ is evident, implying that the ability of the model organism to respond to weak gradients could bethe key to understanding the persistent patterns observed in nature.Turning towards the long term behavior seen in Sections 4.2 and 5.2, by im-plementing only run-and-tumble and autochemotaxis, we gather rich behaviorin one- and two-dimensions that is very similar to the the pattern formationwe see in our model organism (see Fig. 2). From simulations seen in Figs. 5and 9, we observe several aggregations initially form and merge, but as ag-gregations coarsened and expanded, merging events occurred on much largertimescales. Specifically in the two-dimensional case, simulations such as Figs.10a and 10c, which show the rapid organization of aggregations and subse-quent long-time saturation towards larger, non-circular aggregations, reflectthe structure and evolution patterns we see in experiments highlighted in Sec-tion 2. As expected, we are able to predict that the chemical concentration ishighest near the peaks of the aggregations. However, the shape and evolutionof these plankton aggregations and chemical concentrations are dependent onthe properties of the deposition function and δ .We studied three different deposition functions and saw very different re-sults for all three. By comparing deposition functions that admit similar mostunstable wave numbers, we were able to focus on the non-linear developmentafter the initial perturbation becomes apparent. In both the one- and two-dimensional cases, the constant function, f , and the linear switch function, f , are the most comparable to the evolution dynamics seen in our simula-tions when light saturated the environment. The switch function, f , forcesthe plankton to aggregate more slowly and move on a completely differenttimescale than the other deposition functions. Even though f produced moreunique aggregations, this may accurately reflect the speed of merging whenthere is no light as in Fig. 1b.Lastly, our long time simulations varying δ in both the one- and two-dimensional cases stress the significance of a plankton’s chemotactic sensi-tivity. As δ increases, the aggregations are much more spread out and lessdefined. Compare these results to experimental long term behavior as in Fig.2e. The aggregations after 60 minutes are relatively varied, as some are denselypopulated whereas others are spread out much more. This may imply that theplankton’s chemotactic sensitivity is influenced by biological mechanics relatedto a plankton’s response towards light. We conducted several experiments to understand aggregation patterns of phy-toplankton using our model organism
Heterosigma akashiwo . From these ex-periments, it is evident that chemotaxis plays an important role in planktonaggregations, as these aggregations form in absence of external fluid flow andlight. To describe this system, we have constructed a one- and two-dimensionalpartial differential equation model of plankton with run-and-tumble motionand autochemotaxis. We explored the complex dynamics of the models by an-
Run-and-Tumble Model with Autochemotaxis 25 alytically describing the linear stability of constant steady states and numer-ically analyzing the nonlinear saturation and long time behavior. By varyingparameters, we investigated the diverse dynamics of this model from stabilityanalysis to structural differences at steady state. Since the deposition func-tions for plankton are unknown, we chose three simple deposition functionsto simulate. Through simulations, these functions revealed striking differencesin chemical and plankton aggregation patterns, some of which are similar toexperimental observations.There is still much to be explored in this line of research. There are severalsignificant biological and physical mechanisms, such as plankton photosynthe-sis, vertical migration towards and away from light, and a non-uniform turningkernel T ( θ (cid:48) , θ ) which is biologically relevant, that need to be incorporated tounderstand plankton aggregation formation (Chen et al. 2018). Extending thismodel to three dimensions may allow us to understand the extreme change inpattern formation seen in Fig. 1a. We also will conduct more experimentationwith collaborators to find parameters and deposition function which describesthe behavior of Heterosigma akaswhio . There are inherent structural differ-ences in the two-dimensional case for varied deposition functions as seen fromthe graphs of S in Fig. 10f. In future work, we look to recreate S t fromvideos of experiments we have conducted and match the experimental S t withvaried deposition functions and parameters in our numerical simulations. Acknowledgments
We would like to acknowledge Prof. Cathy Coyne and Monica Acerbi for pro-viding us with
Heterosigma akashiwo and food and advice for keeping themalive. We also acknowledge the Center for Computational Biology, directed byDr. Michael Shelley, for their expertise in biological mechanisms underlyingthe plankton. As well, thank you to undergraduates Claire Lubash and DianaLi for their help in laboratory experimentation and numerical simulations.
Declarations
Conflicts of interestThe authors declare that they have no conflict of interest.Code AvailabilityFor the open source code which produces the two-dimensional numerical sim-ulations, see the GitHub repository https://github.com/Louminator/Plankton_Signal_RT . For the code which produced all figures and videos in this paper, see theGitHub repository https://github.com/ricknussell31/AutochemotaxisFiles .For all videos in Online Resources, please see this link https://drive.google.com/drive/folders/1AXWG54FCohTtUPOxIOTmcmddl9CVx9-F?usp=sharing . References
Bearon RN, Pedley TJ (2000) Modelling Run-and-Tumble Chemotaxis in aShear Flow. Bull of Math Biol 62:775–791, DOI https://doi.org/10.1006/bulm.2000.0178Berdalet E, Fleming LE, Gowen R, Davidson K, Hess P, Backer LC, MooreSK, Hoagland P, Enevoldsen H (2016) Marine harmful algal blooms, hu-man health and wellbeing: Challenges and opportunities in the 21st cen-tury. J Mar Biol Assoc UK 96(1):61–91, DOI https://doi.org/10.1017/S0025315415001733,
Chen X, Wu Y, Zeng L (2018) Migration of gyrotactic micro-organisms inwater. Water 10(10):1455, DOI https://doi.org/10.3390/w10101455Fu FX, Zhang Y, Warner ME, Feng Y, Sun J, Hutchins DA (2008) A com-parison of future increased CO2 and temperature effects on sympatric Het-erosigma akashiwo and Prorocentrum minimum. Harmful Algae 7(1):76–90,DOI 10.1016/j.hal.2007.05.006Fu FX, Tatters AO, Hutchins DA (2012) Global change and the future ofharmful algal blooms in the ocean. Mar Ecol Prog Ser 470:207–233, DOIhttps://doi.org/10.3354/meps10047Hallegraeff GM, Anderson D, Cembella A (2004) Manual on Harmful MarineMicroalgae. UNESCO, ParisHara Y, Chihara M (1987) Morphology, ultrastructure and taxonomy of theraphidophycean alga Heterosigma akashiwo. Mar Biol 100(2):151–163, DOIhttps://doi.org/10.1007/BF02488320, URL http://link.springer.com/10.1007/BF02488320
Harvey EL, Menden-Deuer S, Rynearson TA (2015) Persistent intra-specificvariation in genetic and behavioral traits in the raphidophyte, Heterosigmaakashiwo. Front Microb 6(NOV), DOI https://doi.org/10.3389/fmicb.2015.01277Hennige SJ, Coyne KJ, Macintyre H, Liefer J, Warner ME (2013) The Photobi-ology of Heterosigma akashiwo. Photoacclimation, Diurnal Periodicity, andits Ability to Rapidly Exploit Exposure to High Light. J Phycol 49(2):349–360, DOI https://doi.org/10.1111/jpy.12043Hillen T (2002) Hyperbolic Models for Chemosensitive Movement. MathModel Methods Appl Sci 12(07):1007–1034, DOI https://doi.org/10.1142/S0218202502002008
Run-and-Tumble Model with Autochemotaxis 27
Hillen T, Painter KJ (2009) A user’s guide to PDE models forchemotaxis. J Math Biol 58(1-2):183–217, DOI https://doi.org/10.1007/s00285-008-0201-3Kawaguchi S (2011) Chemotaxis-growth under the influence of lateralinhibition in a three-component reaction-diffusion system. Nonlinearity24(4):1011–1031, DOI https://doi.org/10.1088/0951-7715/24/4/002Keller EF, Segel LA (1971a) Model for chemotaxis. J Theor Biol 30(2):225–234, DOI https://doi.org/10.1016/0022-5193(71)90050-6Keller EF, Segel LA (1971b) Traveling bands of chemotactic bacteria: A theo-retical analysis. J Theor Biol 30(2):235–248, DOI https://doi.org/10.1016/0022-5193(71)90051-8Lewitus AJ, Horner RA, Caron DA, Garcia-Mendoza E, Hickey BM, HunterM, Huppert DD, Kudela RM, Langlois GW, Largier JL, Lessard EJ,RaLonde R, Jack Rensel JE, Strutton PG, Trainer VL, Tweddle JF (2012)Harmful algal blooms along the North American west coast region: His-tory, trends, causes, and impacts. Harmful Algae 19:133–159, DOI https://doi.org/10.1016/j.hal.2012.06.009Lushi E, Goldstein RE, Shelley MJ (2012) Collective chemotactic dynamicsin the presence of self-generated fluid flows. Phys Rev E 866318(87), DOIhttps://doi.org/10.1103/PhysRevE.86.040902Riebesell U (2008) Climate change: Acid test for marine biodiversity. Nat454(7200):46–7, DOI https://doi.org/10.1038/454046aSaintillan D (2018) Rheology of Active Fluids. DOI https://doi.org/10.1146/annurev-fluid-010816-060049Seymour JR, Ahmed T, Stocker R (2009) Bacterial chemotaxis towards theextracellular products of the toxic phytoplankton Heterosigma akashiwo. JPlankton Res 31(12):1557–1561, DOI 10.1093/plankt/fbp093Sheng J, Malkiel E, Katz J, Adolf J, Belas R, Place AR (2007) Digital holo-graphic microscopy reveals prey-induced changes in swimming behavior ofpredatory dinoflagellates. Proc Natl Acad Sci USA 104(44):17512–7, DOIhttps://doi.org/10.1073/pnas.0704658104Sheng J, Malkiel E, Katz J, Adolf JE, Place AR (2010) A dinoflagellate ex-ploits toxins to immobilize prey prior to ingestion. Proc Natl Acad Sci USA107(5):2082–7, DOI https://doi.org/10.1073/pnas.0912254107Stocker R, Seymour JR (2012) Ecology and Physics of Bacterial Chemotaxisin the Ocean. Microbiology and Molecular Biology Reviews 76(4):792–812,DOI https://doi.org/10.1128/mmbr.00029-12Visser AW, Thygesen UH (2003) Random motility of plankton: Diffusive andaggregative contributions. J Plankton Res 25(9):1157–1168, DOI https://doi.org/10.1093/plankt/25.9.1157
A Laboratory experiments
The
Heterosigma akashiwo were collected from the Delaware Bay in Lewes, Delaware andkept in beakers filled with seawater. The water was infused with nutrients with a combination8 Nicholas J. Russell and Louis F. Rossiof 1 mL of NaNO3, 1 mL of NaH2PO4 · H2O, 1 mL of Na2SiO3 · . × . × × ×
11 centimeters. For the light source, we used an 80 watt incandescent bulb thatattains approximately 1000 lumens, and we changed the angle and distance of the lightsource for the various experiments. The light source was 7-13 cm away from the dish andwas at a 45 ◦ -70 ◦ angle of incidence. To capture the movement of the plankton, we utilizedan Allied Vision Mako G-30 camera equipped with Vimba Viewer software to take photosevery 5 seconds over the duration of the experiments. Windows Media Player was used tostitch the photos together and construct the movies. For videos, visit the YouTube page forthe University of Delaware Math Plankton Team. B Proofs of Propositions
B.1 Proposition 1
Proof
To analyze the roots of (17) as k → ∞ , we divide (17) by k and define α = k − toobtain α λ + (cid:2) d α + (1 − d ) α (cid:3) λ + (cid:2) ( d + 1) α − d α (cid:3) + (cid:18) d − (cid:18) d + cd δ (cid:19) α (cid:19) = 0 . (40)As α → y = α λ , and substituting this y representation into (40), wederive the cubic equation y + (cid:2) d + (1 − d ) α (cid:3) y + (cid:2) ( d + 1) α − d α (cid:3) + (cid:18) d α − (cid:18) d + cd δ (cid:19) α (cid:19) = 0 . (41)We assume that y is a solution of the form y = ∞ (cid:80) n =0 b n α n , where b n ∈ C . Substituting into(41) and setting the terms corresponding to similar powers of α equal to 0, we obtain α : b + d b = 0 ⇒ b = 0 , , − d α : b (3 b + 2 d b ) = 0 α : b (3 b − d −
1) + b (1 + 3 b + d + 2 d b ) + d (cid:0) b (cid:1) = 0 α : b + 3 b b + 2 d b b + ( d + 1) b + (6 b b + 2 d b − d ) b ) b = 0 . Taking b = 0, the α equation gives no information about b . We then use the α equationto arrive at d (1 + b ) = 0 = ⇒ b = − i, i. (42)Taking b = − i and using the α equation, we get i − id b − id − i = 0 = ⇒ b = − . (43)Now taking b = i and using the α equation, we get − i + 2 id b + 2 i d + i = 0 = ⇒ b = − . (44) Run-and-Tumble Model with Autochemotaxis 29Therefore, the three roots are y = − d + o ( α ) ⇒ λ = − d k + O ( k ) y = − iα + − α + o ( α ) ⇒ λ = − ik −
12 + O (cid:0) k − (cid:1) y = iα + − α + o ( α ) ⇒ λ = ik −
12 + O (cid:0) k − (cid:1) . Since Re( λ ) → −∞ , Re( λ ) → − , and Re( λ ) → − as k → ∞ , then k u → − . (cid:117)(cid:116) B.2 Proposition 2
Proof
Recall that R ( k ) = max { Re ( λ i ) } i =1 . Setting k = 0 in (17), we have λ + (1 − d ) λ − d λ = 0 = ⇒ λ = 0 , − , d (45)Therefore, R (0) = 0 or R (0) = d , depending on the sign of d . In order to find the maximumof R ( k ), we find the k such that ∂λ∂k = 0. Using implicit differentiation on (17), we obtain(46) ∂λ∂k (cid:0) λ + 2 λ ( d k − d + 1) + ( d + 1) k − d (cid:1) + (cid:18) d kλ + 2 k ( d + 1) λ + 4 d k − k (cid:18) d + cd δ (cid:19)(cid:19) = 0 . When ∂λ∂k = 0, we get2 k (cid:18) d λ + ( d + 1) λ + 2 d k − (cid:18) d + cd δ (cid:19)(cid:19) = 0 . (47)Since k ≥
0, (47) implies that k = 0 is a critical point and therefore, R (cid:48) (0) = 0. To classifythis critical point, we seek to find R (cid:48)(cid:48) (0). Solving for the second derivative at k = 0 usingimplicit differentiation of (46), we obtain ∂ λ∂k (cid:12)(cid:12)(cid:12) k =0 = − (cid:16) d + cd δ − d λ (cid:17) λ + 2 λ (1 − d ) − d . (48)If R (cid:48)(cid:48) (0) >
0, then by Proposition 1 and the fact that R ( k ) is a continuous function, theremust be a k m such that R (cid:48) ( k m ) = 0 and R ( k m ) >
0. This will give us the bounds in (19).To this end, we first assume that d >
0. Using λ = d , (48) becomes ∂ λ∂k (cid:12)(cid:12)(cid:12) k =0 ,λ = d = R (cid:48)(cid:48) (0) = 2 (cid:16) d + cd δ − d d (cid:17) d + d . (49)To have R (cid:48)(cid:48) (0) >
0, the parameters must satisfy (cid:34) d − d (cid:32) − (cid:114) d d cδ (cid:33)(cid:35) (cid:34) d − d (cid:32) (cid:114) d d cδ (cid:33)(cid:35) < . (50)Since d >
0, (50) implies d < d (cid:20) (cid:113) cd d δ (cid:21) , which gives the upper bound of(19).Now we assume that d < R (0) = 0. Substituting λ = 0 into (49), we obtain ∂ λ∂k (cid:12)(cid:12)(cid:12) k =0 ,λ =0 = R (cid:48)(cid:48) (0) = 2 (cid:16) d + cd δ (cid:17) − d . (51)This implies that d > − cd δ will make R (cid:48)(cid:48) (0) >
0, which gives the lower bound in (19). (cid:117)(cid:116)
B.3 Proposition 3
Proof
Consider the PDE for the chemical concentration, c , defined in (7a) where ρ ( x , t ) = N p (cid:80) i =1 δ ( x − x i ( t )). We solve this system on a periodic domain D = [0 , (cid:96) ] × [0 , (cid:96) ] and define C ( t ) as the total chemical over time (see (37)). Integrating (7a) over D , noting that f ( c ) = f ( c ) = γ , and using the steady state d c = f ( c ) ρ , we obtain (52a) (cid:90)(cid:90) D c t d x = (cid:90)(cid:90) D [ d ∆ c − d c + f ( c ) ρ ] d x (52b) C (cid:48) = d C + (cid:90)(cid:90) D f ( c ) ρ d x (52c) C (cid:48) − d C = ρ(cid:96) N p N p (cid:88) i =1 f ( c ( x i ( t ))) (52d) C (cid:48) − d C = ρ(cid:96) N p ( γN p ) (52e) C (cid:48) − d C = ργ(cid:96) (52f) C ( t ) = ργ(cid:96) d + A e − d t , (52g) C ( t ) = cγ(cid:96) f ( c ) + A e − d t (52h) C ( t ) = c(cid:96) + A e − d t where A is a constant and ρ(cid:96) N p in (52c) is the plankton density coefficient described in (57).Solving for A using the initial condition of C (0) = C , the solution to the ODE becomes(53) C ( t ) = c(cid:96) + ( C − c(cid:96) ) e − d t . Therefore, as t → ∞ , C ( t ) → c(cid:96) and if C = c(cid:96) , C ( t ) = c(cid:96) for all t . (cid:117)(cid:116) C One Dimensional Numerical Simulations
C.1 Numerical Methods
In this section, we discuss the details of the numerical scheme utilized in solving (11a)-(11b)on the interval [0 , (cid:96) ] with periodic boundary conditions ρ (0 , t ) = ρ ( (cid:96), t ) , ρ x (0 , t ) = ρ x ( (cid:96), t ) , c (0 , t ) = c ( (cid:96), t ) , c x (0 , t ) = c x ( (cid:96), t ) (54)We define D as the C m × C m spectral Chebyshev differentiation matrix on the interval [0 , (cid:96) ]with periodic boundary conditions, where C m ∈ Z + . We also denote ρ n and c n as the n thtime step of the plankton density and chemical concentration, ∆ t as the temporal mesh size,and ∆ x = C m (cid:96) as the spatial mesh size.For (11a), we utilize a pseudospectral method: a Crank-Nicolson method for time-stepping and a Chebyshev spectral method for the spatial mesh. To deal with the non-linearity of autochemotactic, we use a forward Euler method solely for that term. Thescheme can be written as follows: c n +1 − c n ∆ t = d D (cid:2) c n + c n +1 (cid:3) − d (cid:2) c n + c n +1 (cid:3) + f ( c n ) ρ n Run-and-Tumble Model with Autochemotaxis 31(55) c n +1 = (cid:2) I − (∆ t/ (cid:0) d D − d I (cid:1)(cid:3) − (cid:2)(cid:0) I + (∆ t/ (cid:0) d D − d I (cid:1)(cid:1) c n + ∆ tf ( c n ) ρ n (cid:3) , where I is the identity matrix of appropriate size. For (11b), we utilize a pseudospectralmethod: a leap-frog method for the temporal mesh and a Chebyshev spectral method forthe spatial mesh. The scheme can be written as follows: ρ n +1 − ρ n + ρ n − (∆ t ) + ρ n +1 − ρ n − t = D ρ n − D (cid:32) Dc n (cid:112) ( Dc n ) + δ ρ n (cid:33) ρ n +1 = 11 + ∆ t/ (cid:40) ρ n + ρ n − (∆ t/ −
1) + (∆ t ) (cid:34) D ρ n − D (cid:32) Dc n (cid:112) ( Dc n ) + δ ρ n (cid:33)(cid:35)(cid:41) (56)To initialize the system, we define values for d , d , δ, (cid:96), γ, c , T , C m , ∆ t , and f ( c ). We ini-tialize c and ρ as constant solutions which satisfy the constant steady state solution of (11a)-(11b), which is d c = f ( c ) ρ . For simulations in Section 4.2, we utilize (cid:96) = 5, ∆ t = 0 . C m = 301, which is spatially and temporally stable. C.2 Parameters
Table 3 shows all parameters used in the one dimensional simulations in Section 4. Detailsabout a specific variable or notation can be found in Tables 1 and 2.Figure d d (cid:96) δ f ( c ) γ c T c Fig. 4a Var. 1 – 0.01 f f f f f f f f f f Table 3
Parameters used in the 1D analysis and simulations discussed in Section 4. “Var.”means the variable was varied and “–” denotes unused variables.
D Two Dimensional Numerical Simulations
D.1 Methods
In this section, we describe the numerical scheme used to solve the two-dimensional modelfrom the initial discussion in Section 5.2. Unlike the 1D scheme, we discretize the densityfield as a linear combination of moving particles representing an ensemble of plankton.Therefore, we introduce the plankton density coefficient, ρ d . Recall that ρ is the constantsteady state solution to the 2D system, and define the particle density as N p (cid:96) . In order2 Nicholas J. Russell and Louis F. Rossito scale the autochemotactic term appropriately, each of the N p particles must representa small aggregation of plankton. To represent the specified density on the computationaldomain using N p particles, we require that ρρ d = N p (cid:96) = ⇒ ρ d = ρ(cid:96) N p . (57)Thus, for the autochemotactic term, each particle deposits chemical with strength ρ d f ( c ).The chemical field is approximated on a mesh. W we use a periodic [0 , (cid:96) ] × [0 , (cid:96) ] domain,and we define (cid:101) D as the 2 C m × C m two-dimensional spectral Chebyshev differentiationmatrix on the domain [0 , (cid:96) ] × [0 , (cid:96) ] with periodic boundary conditions. We denote ρ n and c n as the n th time step of the plankton density and chemical concentration, where ρ n is asummation of point sources ρ n = N p (cid:88) i =1 δ ( x ni ) , (58)and x ni = (cid:104) x ni , y ni (cid:105) is the position of the i th plankton at the n th time step. We define ∆ t as the temporal mesh size, ∆ x = ∆ y = C m (cid:96) as the spatial mesh size, and x (cid:96) and y (cid:96) as thespatial mesh for the chemical concentration, c , in the x and y directions, respectively.To initialize the system, recall that the constant steady state solution to (7a)-(7c) is d c = f ( c ) ρ . We select a c , d , and f ( c ) and then solve for ρ = d cf ( c ) . All N organisms arerandomly distributed throughout the domain and given a random orientation θ i .For the ( n + 1)th iteration, we complete the following steps:1. Decide if plankton will tumble to a new direction using the probability described in(34). If a plankton does tumble, their new direction is selected randomly from a uniformdistribution θ n +1 i = [0 , π ).2. Move all N p plankton to their new position using (35), which updates ρ n to ρ n +1 .3. Compute S n +1 i := f (cid:16) c n (cid:16) x n +1 i (cid:17)(cid:17) , where x n +1 i is the i th plankton’s position calculatedfrom step 2. We compute c n (cid:16) x n +1 i (cid:17) , the chemical concentration at a particle position,using a bivariate spline interpolation.4. Deposit the chemical. In order to do deposit the chemical, we construct the matrix S n +1 = N p (cid:88) i =1 S n +1 i πσ exp (cid:34) ( x n +1 i − x (cid:96) ) + ( y n +1 i − y (cid:96) ) σ (cid:35) , (59)where the variance is σ = K ∆ td and K is an suitable constant. In the simulations,we have chosen K = 3.5. Evolve the chemical concentration using a Crank-Nicholson method, similar to the one-dimensional case in (55): c n +1 = (cid:20) I − ∆ t (cid:16) d (cid:101) D − d I (cid:17)(cid:21) − (cid:20)(cid:18) I + ∆ t (cid:16) d (cid:101) D − d I (cid:17)(cid:19) c n + ∆ tρ d S n +1 (cid:21) , (60)where I is the identity matrix of appropriate size.For all simulations in Section 5, we utilize ∆ t = 0 .
2, ∆ x = ∆ y = 0 . C m = 200,and (cid:96) = 10. Using these parameters along with all others used in simulations described inSection 5.2, the largest wave number we are able to resolve has a magnitude of | k |≈ D.2 Parameters