MODELLING SPATIO-ANGULAR PATTERNS IN CELL BIOLOGY by ALEXANDER MOGILNER M.Eng., Ural Polytechnical Institute, Russia, 1985 Cand.Sci., Institute for Metal Physics, Russia, 1989 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Department of Mathematics) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September 1995 ©Alexander Mogilner, IH5" In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of / 1 C<-X KcpVL fr3L C C$ The University of British Columbia Vancouver, Canada Date 2.C So^ i ^^ -^ r /33S DE-6 (2/88) Abstract In this thesis I investigate the formation of spatio-angular order in populations of inter-acting individuals. Mathematical models are used to describe distributions of orientation, and how these evolve under different assumptions about the turning behaviour of indi-viduals. Three types of models are introduced: in one, turning is described by a gradual direction shift (rather than abrupt transition). In two other variants, interactions between individuals cause abrupt turning and mutual displacements. Linear stability analysis and weakly non-linear (synergetics) analysis of interacting modes are used to understand the nature and stability properties of the steady state solutions. I investigate how nonhomo-geneous pattern evolves close to the bifurcation point and find that individuals tend to cluster together in one direction of alignment. A special limiting case of nearly complete alignment occurs when the rotational dif-fusion of individual objects becomes very slow. In this case, motion of the objects is essentially deterministic, and individuals tend to gather in clusters at various orienta-tions. To understand analytically the behaviour of the deterministic models, I represent the angular distribution as a number of 8-\ike peaks. For weak but nonzero angular diffusion, the peaks are smoothed out. The analysis of this case leads to the study of a singular perturbation problem. Next, I generalize the models to examine the dynamic behaviour of ensembles of individuals whose interactions depend on angular orientations as well as spatial positions. I show how processes such as mutual alignment, pattern formation, and aggregation can be described by sets of partial differential equations containing convolution terms. Such n models appear to contain a rich diversity of possible behaviour and dynamics, depending on the details of the convolution kernels involved. The analysis of the equations, and predictions in several test cases are presented. Finally, I introduce a partial integro-differential model of aggregation phenomena. The model can be derived from a system of chemotactic partial differential equations (PDE). A linear stability analysis reveals that under certain conditions, periodic patterns may emerge. A weakly non-linear analysis predicts stripe patterns. In the limit of slow diffusion I use a Lagrangian approach to derive and analyse the gradient system and obtain periodic peak-like patterns. The possibility of existence of other types of patterns is discussed. I consider the application of the model to a phenomenon known as rippling in a swarm of social bacteria observed in Myxobacter ia . in Table of Contents Abstract ii Table of Contents iv List of tables vii List of figures viii Acknowledgements x 1 Introduct ion 1 1.1 Aggregation and Alignment Phenomena 1 1.2 Local and long-range interactions 6 1.3 Models for aggregation phenomena 8 1.4 Limitations of pre-existing work 10 1.5 Goals and questions to be addressed 11 2 Mode l s of self-al ignment phenomena in cell biology 16 2.1 Introduction 16 2.2 Description of the Models 18 2.2.1 Model I: Free and bound cells 18 2.2.2 Model II : Objects subject to alignment force 20 2.2.3 Model III: Instantaneous alignment 22 2.3 Summary of the behaviour of Model I 24 iv 2.4 Generalization to three-dimensional rotations 26 2.4.1 Properties of the Laplacian, of convolutions and their kernels in 3D in spherical coordinates 28 2.5 Linear Stability Results in 3D 30 2.6 Bifurcation Analysis of Model I in 2D 34 2.7 Bifurcation Analysis of Model I in 3D 39 2.8 Discussion 44 3 Limit ing case of nearly comple te al ignment 46 3.1 Introduction: the Peak Ansatz 46 3.2 Application of the Peak Ansatz to Model I 49 3.2.1 Reduction of Model I to discrete equations governing the peaks . 51 3.2.2 Stability for the case of two competing peaks 53 3.2.3 Stability for the case of a single peak 55 3.2.4 The case of three interacting peaks 55 3.3 The form of the peaks in Model I (weak angular diffusion) 57 3.4 Application of the Peak Ansatz to Model II 59 3.5 The shape of peaks in Model II (weak angular diffusion) 62 3.6 Application of the Peak Ansatz to Model III 64 3.7 Numerical experiments for Model I 65 3.8 Discussion 67 4 Spatio-angular order in populat ions of self-aligning objects 71 4.1 Introduction 71 4.2 Description of the Models 72 4.2.1 Model I: Free and bound cells 73 v 4.2.2 Model II : Objects subject to alignment and aggregation forces . . 76 4.2.3 Model III: Instantaneous aggregation and alignment 77 4.2.4 Eigenfunctions of the Operators 78 4.3 Linear stability analysis 79 4.4 Discussion 83 5 Partial integro-differential mode l for aggregation in a bacterial swarm 88 5.1 Introduction 88 5.1.1 Swarming behavior of Myxobacteria 88 5.1.2 Rippling 91 5.2 Non-local aggregation model 92 5.2.1 ID model 92 5.2.2 2D Model 95 5.3 Travelling wave solution 99 5.4 ID Linear and non-linear stability analysis 102 5.4.1 Linear stability analysis 102 5.4.2 ID Non-linear analysis 106 5.5 2D Linear and non-linear stability analysis 108 5.5.1 Linear stability analysis 108 5.5.2 2D Non-linear analysis 114 5.6 Discussion 117 6 Discuss ion 119 6.1 Results 119 6.1.1 Significance of the results 121 6.2 Open questions 121 vi 6.3 Some other applications 124 6.4 Some final notes 126 Bibl iography 128 A Bifurcations in an angular space 135 A.l Operators and eigenfunctions in 3D 135 A.2 Stability of the homogeneous steady state in 3D 136 A.3 Non-linear analysis 138 A.3.1 Synergetic analysis of Model I in 2D 138 A.3.2 Bifurcation Analysis of Model I in 3D 140 A.4 Results for the Actin-Binding Kernel 141 B Calculat ions for the Peak Ansatz 144 B.l Stability of two peaks 145 B.2 Stability of one peak 146 B.3 Stability of three peaks 148 C Calculat ions for the spatio-angular case 152 D Supplement to the linear stabil i ty analysis of Chapter 5 157 D.l Properties of the bifurcation function (5.176) 157 D.2 Properties of the LHS of the instability criterion (5.188) 158 vn List of tables 1 Table I. Summary of the results of linear stability analysis for the spatio-angular case 160 2 Table II. Summary of the results of linear stability analysis for the bacterial aggre-gation model 161 Vlll List of figures la-c The kernels K used to represent contact alignment phenomena 162 2a-c The expression Kn(l — Kn) together with the left hand side of the instability cri-terion for Model 1 163 3a-d Plots of the functions present in the instability criterion for Model II...164 4a-f Plots of the functions present in the instability criterion for Model II...165 5a-b The bifurcation diagrams (angular case) 166 6a-c Types of angular patterns in the 2D case 167 7a-c Types of angular patterns in the 3D case 168 8a-b The evolution of one peak (single-humped kernel, Model I) 169 9 The evolution of two peaks (single-humped kernel with a large compact support, Model I) 170 10 The evolution of three peaks (single-humped kernel with a small compact support, Model I) 171 11 The evolution of one peak (single-humped kernel with a small compact support; the diffusion coefficient is large, Model I) 172 12a-b The evolution of two unequal peaks (double-humped kernel: the diffusion coefficient is small, Model I) 173 13 The evolution of two equal peaks (double-humped kernel: the diffusion coefficient is large, Model I) 174 ix 14 The complete bifurcation diagram (spatio-angular case, Models I - III) 175 15 The schemes of the evolving spatio-angular patterns (bifurcation scenario, Models I - I I I ) 176 16 The generalized complete bifurcation diagram (spatio-angular case, Models I -III) 177 17a-c Plots of the functions present in the instability criteria for the spatio-angular case (Models I - III) 178 18a-c The linear growth rates of the unstable modes in the spatio-angular case (Models I - I I I ) 179 19a-c The bifurcation diagrams in the spatio-angular case (Models I - III) 180 20 A diagram of the spatial dependence of interaction of a pair of bacteria 181 21 A travelling wave pattern of the bacterial swarm 182 22 The subdivision of the parameter space of the bacterial aggregation model.. 183 23a-f The qualitative forms of the interaction kernels and bifurcation functions (bacterial aggregation model) 184 24a-d A schematic diagram of periodic patterns of stripes and spots in the bacterial swarm 185 x Acknowledgements First and foremost I would like to express gratitude to my research supervisor, Prof. Leah Edelstein-Keshet, for providing continuous advice, support and guidance. I would like to thank Prof. W. Alt and E. Geigant for discussions and very valuable remarks concerning the contents of Chapters 2,3 and 4. I also would like to thank Profs. G. Bard Ermentrout and S. Gueron for their valuable input in Chapters 3 and 5 of this thesis, respectively. I am grateful to my advisory committee for carefully reading drafts of this thesis. The ideas of this thesis were born in the course of discussions with many researchers in Mathematical Biology, of whom I would especially like to thank Profs. L. A. Segel, G. Oster, Drs. G. Civelekoglu, A. Stevens, B. Pfistner, D. Grunbaum, and K. Ladizhansky. I owe a very special gratitude to my friends without whom graduate school would have been no fun at all. Especially, I am delighted to acknowledge the technical help of Alan Boulton. My work was supported by the University Graduate Fellowship from UBC, an NSERC operating grant and a NATO collaborative research grant to L. Edelstein-Keshet. XI Chapter 1 Introduction 1.1 Aggregat ion and Al ignment P h e n o m e n a Pattern formation is one of the most beautiful and fascinating phenomena in the natural sciences. Biology gives a wide and complex spectrum of spatio-temporal patterns on every level, from the molecular to the meta-populational. Here we will be concerned with special cases of spatial patterns evolving in populations of microorganisms. This type of pattern formation ( usually called aggregation phenomena) has been studied most extensively in Mathematical Ecology. One of the illustrative examples of a highly ordered aggregation pattern is a fish school (Okubo, 1980; Huth and Wissel, 1989). Fish equalize their velocities and directions of motion with one another. An important feature of the interaction between members of a school is its finite-range character: the interaction between a given fish and its neighbours is significant only within some limited spatial domain which contains a number of nearest neighbors. Roughly speaking, the interactions may include tendencies to align with the neighbor, to repel at a short distance, to attract at a longer distance. These interactions result in the formation of schools with characteristic shapes. School characteristics may depend on whether the population density is higher than some critical level. Moreover, not only the shape of a school as a whole may be preserved, but also strong correlations can exist between positions of neighboring fish inside the group. The observed structure of the fish school under some conditions resembles the microscopic 1 Chapter 1. Introduction 2 structure of an amorphous solid. Some smaller organisms, for example insects, may organize themselves into patches without revealing any visible microscopic order. Due to complex and poorly understood responses of individuals to their neighbors' behavior, swarms can remain cohesive for many hours despite considerable randomness in the individual movements inside the swarm which would otherwise dissipate them. The integrity of the group is conserved because, on the periphery, average velocities are directed inward (Okubo, 1980). In this Thesis we focus on a yet more microscopic level, that of bacteria. The co-operative behavior of bacteria is drawing increasing attention (Shapiro and Trubatch, 1991; Budrene and Berg, 1991; Ben-Yakov et al., 1994). Bacteria are often viewed as individual organisms. The limitations of this viewpoint were demonstrated by Shapiro and Trubatch (1991). They investigated a colony of Proteus mirabilis which exhibits characteristics of a multicellular organism. These swimming bacteria are involved in a complex cooperative type of motion. The average velocity of a group depends on the density. A single bacterium does not move but "waits" until some passing group incorpo-rates it; strong short range interaction between bacteria are evident. It has been observed that the more actively moving bacteria are more elongated and aligned. Furthermore, experiments with E.coli described by Shapiro and Trubatch (1991), revealed that just before a swarm starts moving, bacteria turn from a round form into an elongated form and align. Shapiro and Trubatch demonstrated the complexity of the sequence of steps involving cell-cell signalling, cellular differentiation, and multicellular coordination of the proliferation and movement. They suggested that the behavior of a single cell must be regulated in an as yet unknown way to produce the observed pattern of movement. Mod-elling is needed to understand these complex developmental patterns. In particular, it may be possible that a small number of biologically realistic assumptions may lead to Chapter 1. Introduction 3 the essential features of such systems. Remarkable quasi-periodic patterns of "radial arrows", concentric rings, "sunflowers", both standing and propagating, transient, and long lasting were reported under a variety of conditions in colonies of E.coli and Salmonella (Ben-Yakov et al., 1994; Blatt and Eisenbach, 1995; Budrene and Berg, 1991). Microscopic movement is hard to observe in this sort of experiment. Cell proliferation and death probably play an important role in this organization. It seems evident that factors such as chemotactic response, limited food availability, and diffusion of the nutrient affect the colony's morphogenesis. Aggregation phenomena take place on an even more microscopic, molecular level in cell biology. Cell motility in many instances depends on the expansion of the cytoskele-ton, a protein filamentous network at the leading margin of a cell (Bray, 1992). This network in many cases consists largely of short rod-like actin filaments cross-linked by a vast variety of actin-binding proteins. When actin-binding proteins are inactive for some reason, actin filaments disperse rapidly, leading to a homogeneous amorphous structure. However, as a rule, actin-binding proteins have an affinity to actin. By associating and dissociating, these proteins mediate an effective attractive interaction between actin fil-aments, leading to a wide spectrum of actin spatial aggregation patterns: dense clumps of actin in the middle of a cell, dense two-dimensional strands of actin etc. (Bray, 1992). In the past, models for population distributions have tended to dwell exclusively on spatial distribution and temporal dynamics. In nature, however, the relative orientations of individuals may have an important influence on their dynamics and interactions, in-creasing the level of complexity of corresponding biological phenomena. Fish schooling is an important example of this: fish not only aggregate and correlate their velocities with each other, they also align. One can think of many other biological examples in which alignment to a common Chapter 1. Introduction 4 direction or a set of common directions occurs. In a flock of flying birds, individuals moving together as a group orient to the same direction of motion. Other typical exam-ples are herds and other highly structured animal aggregates, where many members of the social group align with each other and move in a common direction. This type of behaviour allows the group to stay together, and it is often adopted as a tactic to avoid predation or improve foraging abilities. On the cellular level, a related in vitro phenomenon is the mutual contact and self-alignment in mammalian cells such as fibroblasts growing on a surface (Elsdale, 1972). The corresponding three-dimensional in vivo process plays an important role in morpho-genesis. These cells form dense patches within which they are partially aligned in a sea of rarefied (low density) randomly oriented cells. In the experiments with E.coli, mentioned above, and with Myxobacteria (which will be described in detail in Chapter 5), elongated bacteria align as the bacterial swarm density increases. On the molecular level, self-alignment in the filamentous networks of the cytoskeleton is often observed and plays an important (and not completely understood) role in cellular motility. In one case, that of filopodia protrusion, this network consists of highly aligned unipolar actin filaments (Bray, 1992). Similar types of angular self-organization is found in a number of subcellular structures such as the contracti le ring (a tight band of actin filaments at the equator of a cell which contracts in the process of cell division) and stress fibers (bundles of actin localized in periodic bands throughout the cell). There are also more sophisticated kinds of angular structures in the cytoskeleton, for example, orthogonal filamentous networks (Civelekoglu and Edelstein-Keshet, 1994). In all cases, we are interested in self-organization (formation of persistent spatial and angular structures due to internal interactions within a population, see Haken (1977)) Chapter 1. Introduction 5 of these populations. Individual random motility and turning of organisms in a homo-geneous environment leads to gradual dispersal and does not create any pattern. If the tendency to disperse is affected by population density, transient patterns may change, but they would not converge to any persistent spatial or angular structure. However, when some aggregation mechanisms are present together with dispersal (a certain balance of these factors is required), inhomogeneities in the population density develop. Qualitatively different mechanisms can lead to pattern formation. Spatial variation of the environment on its own may lead to aggregation under certain conditions. In this situation, each member of the population behaves independently, and the structure evolved is not determined entirely by the dynamics of this population, but rather by inhomogeneities of the environment. Genuine aggregation takes place when organisms respond to each other in such a way that an overall tendency to aggregate prevails over the individual tendency to disperse. The order then is not imposed by external bias or forces, but arises as a natural consequence of the interactions, starting from a chaotic or random state. This type of social interaction makes the dynamics of the population density-dependent. This class of phenomena will be discussed in this Thesis. Similarly, the presence of an environmental gradient, or of some polarized background is not a requirement for the types of self-alignment discussed above, but if it is present, it biases the selection of a direction with which to align. A discussion of self-alignment as a mechanism for enhancing chemotactic ability of social organisms such as schooling fish towards weak, noisy gradients is given in Grunbaum (1994). Locust swarms orient strongly to the direction of the wind. Flocks of migrating birds may use the earth's magnetic field as a directional cue. On the microscopic scale, motile mammalian cells (such as fibroblasts) tend to align strongly with grooves on an artificial substrate. All these cases are not directly addressed in this paper. Instead we investigate the innate Chapter 1. Introduction 6 patterns of organization formed by a group of self-aligning individuals. 1.2 Local and long-range interactions Apart from the necessity of taking the angular distributions into consideration in the above cases, there is another feature of the interaction which poses a challenging mathe-matical problem, namely, the non-local character of the interaction. There are well-known examples where non-local effects are known to be important. These include, first, ag-gregation in an insect swarm (Turchin, 1986). If there is an attraction between species in a swarm, and the perception range of an insect is much greater than the average dis-crete movement step, then local dynamical equations based on density and derivatives become inaccurate. Turchin (1986) has demonstrated that swarming of some insects us-ing auditory, or visual cues to orient towards conspecifics (e.g., Euphydryas males) can be described adequately under an assumption of long-range character of an interaction. Second, neurophysiology provides another example in which distributed interactions are important. One neuron influences another with a strength depending on the distance between the neurons. Following this idea, Ermentrout and Cowan (1979) used theory to find patterns in brain tissue which were based on the idea of non-local inter-neuron interactions. When the individual motility of organisms is understood (a formidable task by it-self), an even more complicated question of understanding and describing inter-specific interactions arises. Globally speaking, one of two types of interaction is usually assumed (Grunbaum and Okubo, 1994; Murray, 1989). The first one can be roughly called "di-rect" interaction. It usually means that some non-chemical signal is transmitted from one individual to another; the velocity or position of the recipient is then adjusted according to this signal. An example is the interactions within a fish school, where communication Chapter 1. Introduction 7 is maintained through various ways including pressure waves and visual signals (Okubo, 1980). In cell biology, a so-called contact inhibition of motion is described in many cases. Not much is known about the nature of this mechanism, but there is enough empirical evidence to suggest that cells are able to change their direction and speed as a result of the direct physical contact between membrane-located specialized proteins. These two examples are different in one important aspect: in the first case the interaction has a long-range character, and in the second case it is short-ranged and contact-like. The second type of interaction is "indirect", through some diffusing material sub-stances, which are emitted by the organisms. This mechanism, called chemotaxis (Keller and Segel, 1971), usually means a drift of one organism up the gradient of a chemical substance secreted by another organism. There are a few well-known examples: some insects aggregate under the influence of pheromones; endothelial cells grow towards tu-mor angiogenesis factor, etc. Evidently, this type of interaction should be classified as a long-range one. The presence of non-local interactions in biological systems stems from (among oth-ers) "direct" interactions, when there is no intermediate signalling substances. In such situations one or another dependence of the interaction on the distance between individ-uals can be derived based on experimental data. In the next section, we will see that when incorporated in models, this dependence gives rise to integral terms. Both local and long-ranged interaction may be repulsive or attractive. Dispersal itself can effectively play the role of repulsion. However, there is convincing evidence that a repulsive component of interaction exists in some bacterial colonies (Blatt and Eisenbach, 1995; Shi and Zusman, 1994). The competition for space and nutrients may lead to interactions that are inherently repulsive. The presence of attractive chemotaxis Chapter 1. Introduction 8 was established in the case of the Myxobacteria swarm (Lauffenburger, 1984). Con-vincing arguments were made in that paper that chemoattractants have fast dynamics. Chemorepellents also were found recently in the Myxobacteria swarm (Shi and Zusman, 1994). 1.3 Mode l s for aggregation phenomena Two main mathematical approaches to the quantitative description of aggregation phe-nomena are the Lagrangian and the Eulerian ones (Grunbaum and Okubo, 1994). In the first, one follows the movement of an individual organism. In the Eulerian case, the spatial density of a population is introduced, and a flux of this population is determined by certain spatio-temporal dynamics. The Eulerian approach is more popular and usu-ally leads to (a set of) non-linear PDEs (partial differential equations) (Murray, 1989; Grunbaum and Okubo, 1994; Holmes et al., 1994; Grinrod, 1991). Insect and bacterial aggregation were modelled extensively with this approach. A distinctive feature of PDEs, their local character, provides the power of this mathematical tool. The problem, though, is that the local interactions are usually more appropriate in physics and chemistry than in biology. Non-local interactions call for integro-PDE models (Murray, 1989; Pfistner and Alt, 1989; Kot, 1992; Levin and Segel, 1985; Lewis, 1994; Turchin, 1986; Mogilner and Edelstein-Keshet, 1995c). Usually the aggregation model equations consist of a lin-ear diffusion term and a non-linear advection term, containing integral expressions for the advection velocity. There are a number of mathematical factors calling for non-local models. It is well-known that a single PDE does not necessarily provide the existence of a stable non-homogeneous solution (Turchin, 1986). A single integro-PDE, on the other hand, allows Chapter 1. Introduction 9 non-trivial solutions to exist (Levin and Segel, 1985). Obviously, there are simplifica-tions in both the analytical and numerical treatment when the number of equations is reduced. There is also a more subtle reason. Some cases of population dispersal are de-scribed satisfactorily by non-linear diffusion equations (Murray, 1989; Grinrod, 1991). It would be natural to model aggregation with some sort of "negative" non-linear diffusion equation (Turchin, 1986). Unfortunately, this mathematical problem is not well-posed (Alt, 1985). What can "cure" this situation is the introduction of short-range integral convolution terms. The clear disadvantage of integro-PDEs is that mathematically they are more difficult to treat. Thus, they are still considered as a "luxury" (Kot, 1992). However, the linear stability analysis of integro-PDEs is no more difficult than that of PDE's (Levin and Segel, 1985; Mogilner and Edelstein-Keshet, 1995c). There are some examples of weakly non-linear analyses of integro-PDEs(Mogilner and Edelstein-Keshet, 1995a). In a number of (somewhat artificial) cases, explicit solutions are known (Grinrod, 1991; Mogilner et al., 1995b). Some rigorous results about existence and uniqueness of the patterns given by integro-PDEs were obtained by Nagai and Mimura (1983), Ermentrout and McLeod (1993). A common approach to integro-PDEs is to assume that the solution changes slowly in space. Then one expands the solution in Taylor series and keeps the first few terms. This technique reduces integro-PDEs to PDEs. There are numerous problems with this treatment. When this method is applied, the dispersion relations at large wave vectors are never a good approximations to the correct dispersion relation (Grunbaum and Okubo, 1994). In extreme cases, bifurcation of the homogeneous state for the full integro-PDE model can be missed. Even if the bifurcation is detected, the neglected terms in the Taylor series are usually comparable to the first few terms. The connection between cells, correlating their orientations, and models containing Chapter 1. Introduction 10 integral equations becomes apparent when one considers that a given cell can meet and interact with a neighboring cell of any other relative orientation. Each such event is associated with some probability that one of the two cells will turn and take on a new orientation. To account for all possible occurrences, one needs to sum up the probability of encounter, weighted by the likelihood of turning over all possible angles of contact. This reasoning leads to the formulation of integral equation models. Similar considerations apply in the binding of macromolecules which must take on the appropriate relative conformations to expose active sites to one another. 1.4 Limitat ions of pre-exist ing work In the existing integro-differential models in Mathematical Ecology the integral terms were treated with the help of linear stability analysis, numerical computations (Levin and Segel, 1985; Turchin, 1986) or reduction to the usual differential operators (Lewis, 1994). The apparatus of functional analysis has been used to prove theorems about the general features of the solutions (Kot, 1992). The angular distribution is usually either ignored or treated separately from the spatial one. When an angular patterning is considered, then a spatial order is ignored. For example, the self-alignment of fibroblasts was investigated in (Edelstein-Keshet and Ermentrout, 1990) under the assumption of the spatial homo-geneity of the cellular population. This assumption is limiting from the biological point of view. Another limitation of that work was the two-dimensional character of the model, excluding the application to the three-dimensional case in vivo. Very general models for complex spatio-angular phenomena were introduced for aggregation and alignment in a bacterial swarm (Alt and Pfistner, 1989) and for branching networks (Ermentrout and Edelstein-Keshet, 1991). In the latter, some simulations and limited treatment of the spatio-angular case were done. The orientations of actin filaments which are bound Chapter 1. Introduction 11 and assembled into a scaffolding structure inside the cell by actin-binding proteins is an example leading to an integro-differential spatio-angular model on the molecular level. This problem concerning angular aspects in a two-dimensional situation has been treated in detail in the work of Civelekoglu and Edelstein-Keshet (1994). The question of which type of spatio-angular order occurs in the cytoskeleton is open. 1.5 Goals and quest ions t o be addressed In this Thesis I shall consider populations of interacting cells or molecules in which individuals tend to aggregate and align. The integral kernels and parameters of the models, in principle, should be derived from experimental biological observations about the details of the motion and interactions of the individuals. However, we do not have enough quantitative data and restrict ourselves to the kernels and parameters postulated which are based on qualitative observations. These models will be treated with the help of a whole arsenal of mathematical tools including linear stability analysis, non-linear bifurcation analysis, and singular perturbation analysis in limiting cases. Analytic results will be supported with numerical experiments: numerical solutions of integro-PDEs and cellular automata experiments. In special cases, the explicit solutions will be obtained in order to support the qualitative analysis. In Chapter 2, I will start with the case of self-alignment in a spatially homogeneous population (if spatial diffusion is fast, then an angular pattern would evolve while the spatial distribution remains homogeneous: the objects undergo alignment but they do not aggregate). The linear stability analysis and synergetic (weakly non-linear bifurcation) analysis of interacting angular modes reveal how non-homogeneous pattern in the angular space evolves close to the bifurcation point at high density or small rotational diffusion rate. The individuals tend to cluster together in one direction. I consider a model for Chapter 1. Introduction 12 angular alignment in 3 dimensions . I start with a model for fast alignment which means that the turning times are much shorter than the characteristic relaxation times for the system. This suggests a turning probability and corresponding integral term in the integro-PDE rather than a detailed description of turning temporal dynamics. This model (Model I) was first applied by Edelstein-Keshet and Ermentrout (1990) to populations of fibroblasts. In addition to this model, I also plan to consider two related models. In Mode l II individuals gradually "drift" (rather than "jump") towards new directions of alignment. The equation of the model describes the convectional drift of the objects in angular space causing alignment. The tendency towards order competes with the dispersal influence of the rotational diffusion. This is a more accurate description if interactions involve actual forces exerted by one unit on another. An example is actin-myosin interactions in the cytoskeleton of the cell. M o d e l III includes density-dependent turning rates. This could account for interac-tions in large populations in which the freedom of movement is restricted due to crowding. In all three models, the nonlinear terms are responsible for aggregation and alignment of the objects, while the diffusion terms describe dispersal of the objects in physical and angular space. We expect, then, that with slow enough diffusion rates or at conditions favouring strong interaction, the degree of alignment will be high. In Chapter 3, I will consider in detail a special limiting situation of nearly complete alignment far from the bifurcation when the rotational diffusion rate is very slow. In this case, stationary angular distributions in the form of narrow peaks will be considered. The stability of these distributions will be analyzed and supported by some singular perturbation analysis and numerical experiments. After obtaining the results on the self-alignment, I will turn to the analysis of more Chapter 1. Introduction 13 complex spatio-angular types of order in Chapter 4. The linear stability analysis will provide the existence of two kinds of bifurcations different from the purely angular one. One of them indicates the onset of the aggregation in an angularly random population. Another one corresponds to some more sophisticated spatio-angular types of order. In these chapters, the models have very general mathematical character, relate qual-itatively to a wide class of biological phenomena, and are not concerned with biological details. The goal is to clarify the qualitative character of spatio-angular patterns which are or will be observed in populations of elongated cells and molecules and to connect these patterns to corresponding interactions. Our goal in Chapter 5 is to provide a non-trivial but analytically tractable mathematical model for bacterial swarming. In this chapter we ignore the alignment processes in a swarm because in a dense swarm, the complete alignment due to tight packing was observed. We concentrate on complex spatial aggregation behavior. This chapter is less abstract than the previous ones. It is concerned with the case of a Myxobacter ia l swarm. Myxobacteria - procaryotic gliding cells - achieve close proximity to each other, creating an aligned coherent swarm. In some cases, the density of the swarm exhibits the stable periodic trains of travelling waves. 1 will a t tempt to explain this phenomenon (called "rippling"). An integro-differential model will be derived and treated with the methods similar to the ones used in Chapters 2 and 4. We will conclude the Thesis with Chapter 6, where the results will be discussed from the mathematical and biological viewpoints. Some other applications of the integro-differential models will be mentioned. Also, some open questions will be posed there. My general goal here is to introduce models for the description of aggregation and alignment phenomena in cell biology based on biologically plausible assumptions. I will be interested in cases where a non-local character of interactions is essential. The desired Chapter 1. Introduction 14 result is to obtain examples of biologically meaningful pattern formation in physical and angular space. More specifically, I will pursue the following goals and attempt to answer the following questions: 1. Introduce a set of models for the behavior of randomly turning cells or molecules subject to a density- and angle-dependent non-local interaction. A. What are the conditions for instability of the homogeneous angular distribu-tion? B. What kind of pattern evolves close to the bifurcation point? C. How should one describe angular distribution in a limiting case of nearly com-plete alignment? D. Are these non-homogeneous angular distributions stable? E. Are these results robust? 2. Generalize the angular models on spatio-angular situations of spatially inhomoge-neous populations with non-local angularly and long-ranged spatially interactions. A. What are the conditions for instability of the homogeneous distribution in this case? B. What are possible spatio-angular patterns close to the bifurcations? Which features of spatio-angular interactions are responsible for one or another pat-tern? C. Are these results robust? 3. Derive a non-local aggregation model for a bacterial swarm from the assumption of chemotactic interactions. Chapter 1. Introduction 15 A. How does one explain the long-range travelling wave shape of the swarm? B. What are the conditions for instability of the homogeneous density distribution in the swarm? C. What are the specific features of the interactions responsible for the peculiar short-range "rippling" pattern? The main goal of the integro-PDE approach to aggregation and alignment phenomena is simply to see which set or sets of biologically plausible elementary rules of individual dispersal and interactions can produce the observed patterns. An interesting but as yet unsolved problem is the "inverse" problem: how to predict the form of the interactions from the observed patterns. The biological applications in areas of cell biology such as morphogenesis or wound healing are possible. Chapter 2 Mode l s of self-alignment phenomena in cell biology 2.1 Introduct ion In Chapters 2 and 3 pattern formation phenomena in which patterns of angular (rather than spatial) distributions form are considered. Neglecting the spatial distributions will be justified in Chapter 4. Briefly, if the spatial diffusion rate is effectively greater than the rotational diffusion rate, then a population is well mixed and homogeneous in physical space, while the angular distributions can be non-homogeneous. In this chapter, I present results and analysis of models for contact-induced turning responses and alignment in populations of interacting individuals. Such models describe distributions of orientation, and how these evolve under different assumptions about the turning behaviour of indi-viduals. One of these models was first used to describe interactions between mammalian cells and fibroblasts in Edelstein-Keshet and Ermentrout, 1990 (henceforth abbreviated EKE (1990)). The original model (first introduced in EKE (1990) and called M o d e l I) is generalized to encompass motion in both 2 and 3 dimensions. I broaden the discussion of EKE (1990) to a more general exposition of the problem of orientational order and its mathematical analysis. Two modifications of this model are introduced. In Mode l II individuals gradually "drift" (rather than "jump") towards new directions of align-ment. This is a more accurate description if interactions involve actual forces exerted by one unit on another. An example is actin-myosin interactions in the cytoskeleton of the cell. The third model (Model III), includes density-dependent turning rates. This 16 Chapter 2. Models of self-alignment phenomena in cell biology 17 could account for interactions in large populations in which the freedom of movement is restricted due to crowding. Using linear stability analysis and weakly non-linear analysis of interacting modes I describe the nature and stability properties of the steady state solutions. The evolution of nonhomogeneous pattern close to the bifurcation point is investigated. It will be found that individuals tend to cluster together in one direction of alignment. Model I was first applied in EKE (1990) to populations of fibroblasts, mammalian cells that move by crawling. The cells are polar, having a "front" and a "rear". In the living body, cells navigate through a complex three-dimensional matrix. In artificial culture, the cells adhere to and move along the surfaces of their culture flask (on a 2D surface). The cells interact with one another as follows: If direct contact occurs between a given cell and its neighbor, the given cell will either retract or turn. The neighbor sometimes turns, also. This turning response forms the main phenomenon of interest here. The response depends on the relative orientations of the participating cells (Elsdale, 1972; Elsdale and Wasoff, 1976). The phenomenon at the population level is that cells initially randomly oriented can form patches of alignment called parallel arrays (Elsdale, 1972; Elsdale and Wasoff, 1976). Qualitatively similar (and biologically more significant) phenomena take place in a three-dimensional situation in vivo in the course of certain stages of morphogenesis. Model I consists of a pair of integro-partial differential equations which describe an-gular distributions of cells. Here the assumption that the cells are distributed uniformly with respect to position will be made. Formation of preferred directions of alignment in a 2D region can be described as formation of pattern on a one-dimensional domain with pe-riodic boundary conditions, or more simply, pattern formation on a unit circle. My goal in this chapter is to generalize the analysis in EKE (1990) to angular distributions Chapter 2. Models of self-alignment phenomena in cell biology 18 in 3D space, that is, to pattern formation on the unit sphere. I also investigate the influence of the drift term and the density-dependent term by comparing the first model with Models II and III. These models are suitable not only for cellular phenomena, but for some molecular phenomena. In Civelekoglu and Edelstein-Keshet (1994), structures formed by actin fibers crosslinked by actin-binding proteins were described. This problem also fits into the general scope of formation of orientational order. In the next chapter, we consider a limiting situation of a slow rotational diffusion leading to almost complete alignment. Here we concentrate on a parameter region close to the bifurcation value, so that the amplitude of the patterns are small. 2.2 Descr ipt ion of the Mode l s 2.2.1 Mode l I: Free and bound cells As in EKE (1990), I restrict attention to cells distributed uniformly over space. Two groups of cells are considered: (1) those that are associated in groups ("bound cells"), share a common axis of orientation, and are restricted to motion along a single direc-tion; and (2) those that are isolated ("free cells") and can continually undergo random reorientation. If the medium is isotropic, it is reasonable to assume that the probability of alignment of two cells depends only on the relative angles between the cells. In Model I, it is assumed that changes in orientation take place "instantaneously" when two cells align. The following variables are used in the analysis: Definit ions t = time, 0 = angle of orientation relative to some arbitrary direction (e.g. the x axis), Chapter 2. Models of self-alignment phenomena in cell biology 19 C(6, t) = density of free cells oriented with angle 6 at time t, P(Q,t) — density of bound cells oriented with angle 6 at t ime t, (3 = rate of alignment due to contact, 7 = rate of exchange between bound and free cells, /J, = rate of random turning by free cells. K(0, 9') — K{6 — 6') = angle dependent part of the rate of contact-induced turning from 6' to 0, It is understood that angles are taken in the range —7r < 6 < TT and that functions of 6 are periodic with period 27T. As in EKE (1990) we consider the set of equations f ( 0 , t ) = frCK* C + foPK*C-jP, f£(0, t) = [i/^eC - pxCK * C - p2CK *P + jP, (2.1) where the usual notation for convolution have been used, K * C{9) = j " K{0 - 6')C{6\ t)de'. (2.2) Here A# is the Laplacian operator. In the 2D case, A# = J ^ . In equations (2.1), fj, represents the rate of random turning of free cells and 7 is the rate of shedding of cells from the bound fraction. The convolution term (2.2) captures the following elementary process of alignment: two free cells meet with initial directions 6,6'. With some prob-ability, say a, the cells continue moving with no interaction and directions unchanged. We make an assumption that with equal probability (1 — « ) / 2 per unit time, the cells align to either direction 6 or 6'. (The factor (1 — a ) / 2 is absorbed into the constant /?i.) This assumption is justified by some experimental observations (Elsdale, 1972), but other reasonable assumptions (as, for example, both cells move to (6 + 6')/2) lead to qualitatively similar results. After aligning, the cells stick, becoming part of the class of Chapter 2. Models of self-alignment phenomena in cell biology 20 bound cells. The term fixCK * C (which stands for (3iC{K * C)) gives the rate at which free cells are converted to cells bound at angle 9 through contact with other free cells. The term /?2 PK * C is the rate at which free cells accumulate at angle 9 due to contact with bound cells at that angle. By comparison, the term —ft2CK * P is the rate of loss of free cells from angle 9 due to attraction to bound cells at any other angle. Note that these terms are not symmetric. The term K * C can be viewed as a dimensionless expression which depends on the entire distribution C(.) , but which acts at a given angle 9. We will refer to this expression as the influence of the free cell distr ibution on the direct ion 9. Similarly, we will refer to K * P as the influence of the bound cell d is tr ibut ion on the direction 9. 2.2.2 M o d e l II : Objects subject to a l ignment force In Model II, turning is described as a gradual rotation, rather than a spontaneous jump to a new orientation. This might be appropriate if alignment were mediated by an applied "force" (e.g., elastic tension). I describe the turning motion of the object as a rotation about its center of mass with some angular velocity. The angular velocity will result from a force acting on the individual. In physics, many forces (such as those in electrostatics, but in particular any conservative force) can be expressed in terms of some underlying scalar potential function. I illustrate below how, in Model II, the kernels of the convolutions that describe interactions are most naturally interpreted as such potentials. The angular velocity of an individual is defined to be w = d9/dt. Since I am here describing the model in 2D, to = d9/dt is a scalar. Later the model will be generalized to 3D, in which case to = dQ/dt is a vector, ft would be the angle between two points fi l 5 fi2 on the sphere, and the direction of the vector fi would be tangent to the arc joining Q,i Chapter 2. Models of self-alignment phenomena in cell biology 21 and 02-I assume that the angular velocity is proportional to the corresponding force: w = Fe, (2.3) The proportionality assumption is characteristic of movement through fluid at low Reynolds numbers. (Without loss of generality I take the constant to be 1.) The force Fg is represented as a gradient of some corresponding potential field, W: Fe = VeW, (2.4) The potential at a given angle 6 is created by the sum total of interaction with individuals at all other angles. I assume linear superposition of forces. Then the potential W can be written in the form: W = W*C = fw(0-0')C(0',t)d9'. (2.5) Since each angle is associated with an angular velocity, objects tend to "drift" collec-tively. This drift can be described by a convection term in a canonical way: Je = Ceo (2.6) Including the term cAgC for angular diffusion, leads to a balance equation for the angular distribution as follows: BC — = eAeC-V0-39. (2.7) Putt ing together equations (2.3-6), I arrive at Mode l II BC — = eAeC - Ve • (CVe(W * C)) . (2.8) This equation describes the convectional drift of the objects in angular space towards the points of the highest concentration, causing alignment. This tendency towards order competes with the dispersal influence of the diffusion. Chapter 2. Models of self-alignment phenomena in cell biology 22 2.2.3 M o d e l I I I : I n s t a n t a n e o u s a l i g n m e n t The third model can be viewed as a rough simplification of the second model, if we again drop the local nature of the interactions and introduce instantaneous jumps from an initial orientation to one resulting from an interaction. I use the same notation, C for density, and omit the distinction between free and bound objects, as in the second model. I introduce the following equation: M o d e l I I I : rtC — = C(Q(C)*C) + eAeC. (2.9) Here the integral term is defined as: Q(C) *C= I dO' L{C{6) - C(9'))G{6 - e')C{6', t). (2.10) The kernel is now a product of two functions. The density dependent function L(C(6) — C(8')) reflects the tendency for a bigger cluster to grow at the expense of the smaller cluster. (With slight abuse of the terminology here by "cluster" I mean a group of objects of the same orientation, not close to each other spatially.) The function G describes the angle dependence and has the same meaning and form as K in Model I. This is a rough approximation of the process of fast rotation of a small cluster of objects towards a more slowly moving big cluster and their final merging. I assume that the function L is odd L(—C) = —L(C), monotonic and bounded and L'(0) > 0. That L is odd means that the nonlinear term in equation (2.9) describes either growth or decay of the cluster at 9, depending on relative size of other clusters. Further, symmetry of the function L provides conservation of the total mass of the system. Models I-III share several limitations: Chapter 2. Models of self-alignment phenomena in cell biology 23 1. These models do not describe extremely low densities, for which stochastic processes of cell movement cannot be successfully approximated by continuous PDEs. 2. The models are inappropriate for extremely high densities when topological pack-ing constraints dominate all other effects (see Elsdale and Wasoff (1976); Onsager (1949)). 3. The binding rates of two free versus free and bound individuals are probably dif-ferent, but I have neglected this distinction, taking /?i = /?2, and using the same kernel to describe both events. 4. One of the approximations made is that I do not distinguish between cell clusters of different sizes. To have a detailed account of these, one would have to introduce the functions P2, P3, P4, • • • to denote clusters composed of two, three, four, ...and n cells. The result would be a system of infinitely many equations for C, P2, P3, •.. which would lead to a complicated mathematical problem. Instead of doing so, I define the simplified variable P = J2%>=2 nPn-The detailed behaviour of the models depends in an interesting way on the nature of the kernels K, W, G. These kernels would have distinct properties in each model problem considered, as they are based on details of the biology that are either observed through experiments or conjectured from some knowledge of the system. In the case of fibroblasts, interactions causing alignment are known to be weakest if the cells meet at 90°. (See Figure la) . In the case of actin fibers, crosslinking proteins of various sorts allow fibers to interact and bind at different configurations, including parallel and orthogonal structures (Civelekoglu and EK, 1994). The kernel is different in that case (see Figure lb-c). In parallel interactions, we must still consider a further distinction, namely whether alignment occurs only "head-to-head" or also "head-to-tail". The first Chapter 2. Models of self-alignment phenomena in cell biology 24 case leads to a kernel with a single hump in the domain — IT < 9 < IT (See Figure la) . The second case results in a kernel with a double-hump (See Figure lb) . In this chapter, I will focus attention mainly on the case of the single-humped kernel discussed in EKE (1990), i.e., on the case that cells align only head-to head, and not head-to tail. The case of the double humped kernel is a simple generalization of this situation which is obtained in 2D if angular space is changed from [—7r,7r] to [—|, | ] , and in 3D by an analogous change of angular space, fL Single- and double-humped kernels lead to parallel alignment of cells. In the case of actin, where interactions are strongest between orthogonal fibers, two mutually orthogonal axes of orientation can be formed. The total mass of cells is equally distributed between these two axes. This case can be treated within the framework of a model essentially identical to Model I, provided the form of the kernel is suitably adjusted, as described in Appendix A.4. 2.3 S u m m a r y of the behaviour of Model I We first concentrate on the first model and briefly review results obtained in (EKE, 1990; Civelekoglu and EK, 1994). We make several observations about the equations of Model I: 1. The total mass of bound and free cells at all orientations, M, is a constant. This can be seen by integrating both sides of equations (2.1) and adding the two equations. Thus, the equations admit a conserved quantity. 2. The equations have a homogeneous ^-independent steady state, P, C in which (in the case (3\ = /?2 = /?) P/C = (3M/j. (2.11) Chapter 2. Models of self-alignment phenomena in cell biology 25 3. The equations can be linearized about this steady state (C = C + c, P = P -f p) yielding § f ( M ) = PiCK * c + (32PK * c + C(f3lC + i32P) - 7P, f ? ( M ) = A # - PiCK * c- foCK * p - ( ^ C + fcP)c + 1P. 4. These linearized equations contain two linear operators, (a) the Laplacian in ID with periodic boundaries, and (b) the convolution on the same region. The set of eigenfunctions of the Laplacian in ID, are simply the functions ^n(0) = eine n = 0,L.. . (2.13) or equivalently, sines and cosines of (nQ). The integer n is called the m o d e num-ber. The eigenvalue corresponding to the eigenfunction ipn = emS is — n2. 5. A fundamental fact in this model is that the eigenfunctions of the Laplacian are also eigenfunctions of the linearized integral operator, with the property that K*ipn = Knipn, n-0,1,... (2-14) where K is the Fourier transform of the function K{6). The fact that the Laplacian and the integral operator share a set of eigenfunctions is essential in the analysis of the problem. (If this is not the case, numerical analysis is needed.) 6. The linear analysis of this problem was explained in EKE (1990). Briefly, it was found that the homogeneous steady state could be destabilized by a perturbation for modes n that satisfied the dispersion relation An2 < Kn(l - Kn), (2.15) where A^JM?- <2 '16> Chapter 2. Models of self-alignment phenomena in cell biology 26 In the case of a single-humped kernel, it is found that the wavenumber n = 1 is the first destabilizing mode, whereas if the kernel has two humps, the first unstable mode is n = 2. 7. The case of the actin alignment kernel (Figure lc) is discussed in Appendix A.4. The wavenumber n = 4 is the first destabilizing mode, both in the 2D and the 3D cases. The mathematical treatment of this case is completely analogous to the single and double humped kernels, and therefore we omit detailed development. 2.4 General izat ion to three-dimensional rotations When rotational motion takes place in three dimensions, one must describe the set of possible directions by two angular variables, for example by the angles <j> and 0 used in spherical coordinates. A direction in 3D can be represented by a unit vector. Further, this vector can be identified with a point on the unit sphere. Thus if we were interested in the distribution of directions in a discrete population (composed of, say, M individuals), we could represent this by a distribution of M points on a unit sphere. However, as our models are concerned with continuous angular distributions in 3D, I shall deal with functions on the unit sphere. Thus, alignment of the cells is equivalent to formation of pattern on the unit sphere. (See, for example, Hunding (1982)). I will interchangeably refer to "cells distributed at various orientations" and "cells distributed on the unit sphere". The angular coordinates Q, = (</>,#) will be used to denote position on a unit sphere. These coordinates are defined with respect to a cartesian coordinate system by the usual spherical coordinate transformation (x,y,z) = (cos 0s'm<f), sin#sin</>, cos</>). (2-17) Chapter 2. Models of self-alignment phenomena in cell biology 27 The previous problem can be generalized to 3D geometry by converting the operators to a fully 3D form. As in the 3D case, rotational diffusion can be represented as a random walk in the angular space, i.e., on the surface of a unit sphere. Thus, this process can be described by the (angular part of the) Laplacian operator (Priestly et al., 1975). From here on, when discussing "the Laplacian" I shall refer only to the surface spherical part of the operator. We find that the equations can be written in the form dp at (Q, t) = PiCK *C + 02PK *C--yP, § ( f t , t) = fiAnC - p^CK * C - p2CK *P + jP. (2.18) where now P(fl , t) and C(tt, t) are functions defined on a unit sphere. The Laplacian operator An in surface spherical coordinates is shown in Appendix A. l . The convolutions are now done in 3D spherical geometry, so that K*C=I K(n, n')c(n')dn'. (2.19) (See Appendix A.l for the full <j>, 9 form.) This expression can be interpreted as a clearcut analogue to the convolution in 2D. That is, we can think of K * C as t h e influence of t h e free cell distribution on the direction Q, = (c/>,0). A similar interpretation can be made for K * P. I assume, as in EKE (1990), that P\ — p2. The equations can then be brought into the following form: ^(n,t) = CK*C + PK*C -aP, (2.20) (o i\ — e \r. — r. K ± r.- r K ± p J. nP dt dC-(n, t) = eAC -CK*C-CK*P + aP. where a = j/P, e = a./P (and t —> pt). A similar generalization can be made in the case of each of Models II and III. The equations of these models are easily rewritten in terms of the spherical angle 0 . Chapter 2. Models of self-alignment phenomena in cell biology 28 Linear stability analysis of the model (2.20) leads to a linear integro-partial differential equation problem, which together with the boundary conditions and geometry forms an eigenvalue problem. I consider separately the eigenfunctions and eigenvalues of the two operators that appear in the equations. 2.4.1 Propert ies of the Laplacian, of convolutions and their kernels in 3 D in spherical coordinates The Laplacian operator in the surface spherical geometry has as its eigenfunctions the surface spherical harmonics (SSH) Yn{<j),0) of degree n with AYn = -n(n + \)Yn (2.21) i.e. the corresponding eigenvalues are A„ = — n{n + 1), n = 0 , 1 , 2 , . . . . The form of Yn is given in Appendix A.l and involves a linear combinations of P°, Legendre polynomials of degree n, and P™, associated Legendre functions of degree n and order m. (Macrobert, 1967; Kraut, 1979). Due to symmetries of the sphere, the nth eigenvalue of the Laplacian has a (2n + 1)-fold degeneracy. This fact plays a very important role in the bifurcation analysis of the 3D case. Note that due to this degeneracy of the spherical eigenfunctions, the eigenvalues depend only on n and not on m. (See Macrobert (1967)). The first few Legendre polynomials for the case m = 0, namely, P° , (usually written simply as Pn) are given in Appendix A. l . As in the 2D case, if we assume that the environment is isotropic, the interactions between two cells depend only on the relative angle of contact of these cells. In 3D, the convolution K * C takes the form K*C= I K(n - tt')C(tt')dn', (2.22) Chapter 2. Models of self-alignment phenomena in cell biology 29 where the integral is taken over the surface of the unit sphere S. We can write this convolution in terms of the angle between the two interacting individuals, 7, but it is more convenient to express it in terms of the cosine of this angle, cos 7. Defining r\ = cos 7, and rewriting the kernel as a function of this argument but using the same symbol K for this new function, we have K*C= / K(rj)C((f>f,0')d(cos(f>')dd'. (2.23) We observe that a vector dot product of the two vectors representing the interacting individuals produces a formula for cos 7: Representing the directions of the cells by the two unit vectors n = (cos 6 sin 6, sin 6 sin 6, cos 6), K ' (2.24) n ' = (cos 6' sin cj>', sin 6' sin <f>', cos cf>'). Forming the dot product, and rearranging algebraically leads to the result cos 7 = n • n ' = cos (f> cos (f)' + s in^s in0 'cos (^ — 0'). (2.25) For stability analysis, the fact that we will later use is that the surface spherical harmonics (SSH) are also eigenfunctions of the convolution operator K with the integral kernel K(£l — £l'). This result is analogous to the 2D case. The appropriate eigenfunctions are now the set of functions Pn{r})i n = 0,1,.... This set is a complete orthonormal set of functions, and we can thus write an expansion of the function K(n) in terms of these functions, 0 0 Kir,) = E KPM, (2.26) where K = 2-^± £ K{r,) P0M dr,. (2.27) Chapter 2. Models of self-alignment phenomena in cell biology 30 The inner products of SSH with the Legendre polynomials are also SSH (see Appendix A. l ) . In 3D, we will use the notation Kn = K{n) = 2TT J^ K(rj) P0^) dr,. (2.28) It is evident that Yn((f), 9) are the eigenfunctions of the integral operator K and Kn are the corresponding eigenvalues. It is worth pointing out the similarity between these Kn for a given integer mode number n, in the 3D case and the Fourier transform K(n) for a given integer mode number n in the 2D case. Again we stress the similarity to the 2D case where eigenfunctions of both the Laplacian A# and of the integral operator are the same. We also comment that as in the 2D case, the expressions for transforms of the kernels K(n) are positive for n = l in the single hump case and for n=2 in the double humped case, i.e. 0 < Ka(l),Kd(2) < 1. This follows from the fact that ||P°(a;)|| < 1 which permits integral estimates of the transforms to be made, and from the normalization of the kernels. 2.5 Linear Stabil i ty Resul ts in 3 D Equations (2.20) have a homogenous steady state solution in which the distributions C, P are constant for all directions in 3D. In an analogous way to the 2D case, we consider a perturbation ~p(n,ty c(n,t)_ = ~P~ c_ + ~Po~ Co. K ( 0 ) e At (2.29) Details of the stability calculation are shown in Appendix A.2, and follow closely those in the 2D case. We find that in 3D, instability of the homogeneous distribution occurs at any harmonic n for which the following inequality is satisfied: ^ / 7 An(n + l)<Kn(l-Kn), A = ^ { ^ f . (2.30) Chapter 2. Models of self-alignment phenomena in cell biology 31 Note that this expression is analogous to inequality (2.15) for the 2D case. The coefficient A is identical, but the dependence on the degree of the harmonic n is slightly different. The conclusions from (2.30) are: 1. To get some feeling for this inequality, we chose several representative functional forms for the kernel K(j). In each of these cases, the kernel represents alignment only when objects meet at "small enough" contact angles (angle 7 for which cos 7 > 3/4). Further, the kernels are all single humped, which means that anti-parallel interactions do not cause alignment. We have numerically plotted the functions on both sides of the inequality (2.30) for two values of the constant, A = 0.1 and A = 0.03, see Figure 2, and for kernels of the form f / ( 7 ) , c o s 7 > 3 / 4 ; K = { (2.31) [ 0 , c o s 7 < 3 / 4 . with various functions / , see Appendix A.2. Each of these forms has finite support on a subinterval for which |—yj < | , but the kernels Ki,Ku, Km have different types of discontinuities and shapes. It was shown in EKE (1990) that the number of an unstable mode does not depend on the exact functional form of K{6) in 2D, given its symmetry and type of support. Figure 2a-c show the plots of the expression Kn(l — Kn) as a function of n for three kernels of the form (A.209) and the quadratic curves corresponding to the right-hand side of the inequality and the way that they intersect. 2. From these figures it is clear that as we decrease A, symmetry breaks first as a result of the growth of the first harmonic. If A is much smaller, the second harmonic becomes unstable, then the third, fourth, etc., in succession. However, it is not necessarily true that the mode n = 1 is the one that always breaks stability. Chapter 2. Models of self-alignment phenomena in cell biology 32 In particular, if a kernel contains no components of the eigenfunction corresponding to this mode, one of the higher modes will cause instability. 3. A double -humped kernel occurs if interactions occur at angles both close to zero and close to 180°. The correspondence between the single-humped and the double-humped kernel, Kd, is as follows: { 0, n odd (2.32) A s (n) , n even. This holds provided the shape of Kd is identical to that of Ks but the periodicity is doubled. The result is that odd harmonics can never cause instability in the case of Kd, so that the break of symmetry is caused by the second harmonic. 4. Since A is positive, instability is most likely for low values of the integer n, for example for n = 1. Further, as we show in Appendix A.2, high harmonics cannot destroy stability. 5. In comparing the conditions for instability in 2D and 3D (i.e., the inequalities (2.15) and (2.30)) we observe that in 2D the minimum value of the LHS of the instability criterion (2.15) when n = 1 is mm An2 = A, whereas in 3D the minimum value of the LHS of (2.30) would be min An(n + 1) = 2A. The right hand sides of both inequalities are the same, and in both cases the expression attains a maximum value of max(A'n(l - Kn)) = {. Thus, if A > 1/4 in 2D or A > 1/8 in 3D, the homogeneous state is stable. This shows that the homogeneous state tends to be more stable in 3D. This makes physical sense also because of the larger number of diffusional degrees of freedom in 3D that can destroy order. Stabi l i ty for Mode l s II and III Chapter 2. Models of self-alignment phenomena in cell biology 33 Let us now consider the stability for Models II and III. Both these models contain a single equation, so the analysis is quite simple. In both models the homogeneous solution in which the distribution C is constant for all directions is stationary. We substitute the perturbed distribution C(ft, t) = C + C0zn{Q)eXt (2.33) into the equations (2.8) and keep terms linear in Co- Then instability to growth of the n ' th mode occurs whenever A > 0. Mode l II The instability criterion has the form (both in 2D and 3D): e < CWn (2.34) In Figure 3 this inequality is illustrated graphically for the following step kernels: In 2D: ~, —a < 6 < a W(0) = {2a _ (2.35) 0, otherwise (A2, c o s 7 > | W(Q) = (2.36) (0, cos 7<f Mode l III Using the fact that at small C we have L(C) = r\C + 0(C2), we obtain the following instability criterion: f n2, in 2D vC2(l-Gn)>e\ (2.37) [n(ra + l ) , in 3D The inequality is illustrated in Figure 4 for the same kernels as the ones used for Model II. The statements 1-4 for the case of Model I hold also for Models II and III. In particular, as e is decreasing, stability is broken first by the first harmonic (n = 1) for single humped kernel and by the second harmonic for the double humped one. and in 3D: Chapter 2. Models of self-alignment phenomena in cell biology 34 2.6 Bifurcation Analys is of Mode l I in 2D In order to understand the pattern that arises after stability of the homogeneous distri-bution is lost, we apply an analysis based on the approach outlined by Haken (1977), Friedrich and Haken (1989). The basic idea of this approach is that close to bifurcation, the fastest growing mode essentially controls the amplitudes of other modes which are just becoming unstable. This is the so-called adiabatic or "slaving" principle of synergetics. The rationale for this assumption is as follows: When the amplitudes of all harmonics are very small, the unstable mode grows exponentially. Further, the relaxation time of the stable modes is small relative to the time scale of variation of the unstable mode. This means that the stable modes will quickly reach meta-stable states and thereafter slowly follow the dynamics of the unstable mode. These dynamics lead to an increase in the amplitude of the first mode at the expense of the other modes. We begin with the 2D case having 0(2) rotational symmetry in which the first har-monic is responsible for the stability break. The growth rates of other modes are then neg-ative, (An < 0, n = 2, 3 , . . . ) , and the first mode is just beginning to grow (Ai > 0, \\ ~ 0.). Then this unstable harmonic will be the leading mode, but its linear growth rate magni-tude is small compared to the magnitudes of the decay rates of the other modes close to the onset of instability (|Ai| <C |A,-|, i = 2, 3 , . . . ) . Details of the calculations are given in Appendix A.3. As a first step, we express the state in terms of a superposition of the steady state, the growing first harmonic, and the other harmonics, \P(6,t)) \P) \Vl(t)J .v^iUC)/ In this expansion, the values (£;(£), i]i(t)) are time-dependent amplitudes of the given harmonics, and these depend on the amplitude of the first mode. Note, that £o(0 is the Chapter 2. Models of self-alignment phenomena in cell biology 35 amplitude of the constant component of the free cell angular distribution. It is not zero automatically. This assumption is valid close to criticality where the behaviour of the system is entirely determined by the single unstable mode. The amplitudes of all the other modes are small, so it is possible to form a valid asymptotic expansion for £•(£) and T)i(t) in powers of the unstable mode amplitudes £1, and 771, as shown in Appendix A.3. This approach allows us to eliminate the amplitudes of the stable modes and leads to a closed set of ordinary differential equations for the unstable amplitudes, the so called general ized Ginsburg-Landau equations: ~( j = A M ) +p{tuVi,ti({i,m)MtuVi),} *' = 0 , 2 , . . . , (2.39) where P is a non-linear functional. If we were to analyze the steady state solutions of this set of ordinary differential equations we would get a detailed description of pattern created close to bifurcation. We will look only for the time-independent solutions, and not for the full time behaviour, for which the calculations are forbidding. That is, we set This leads to the equations C(K*C) + P(K*C)-aP = 0, cAeC - C(K * C) - C(K * P) + aP = 0, (where now, in the 2D case, A# = d2/d92). From equations (2.41) we deduce that _C(K*C)_ F~a-{K*C)y ( ' eA8C + V(C) • C = 0, (2.43) { ' a-(K*C) Ka-{K*Cy K ' Chapter 2. Models of self-alignment phenomena in cell biology 36 We now look for a solution of equations (2.43, 2.44) in which the deviation away from the homogenous state is expressed as a superposition of the harmonics. That is, we let C(e) = C + £(0), \t(0)\<C (2-45) where £(#) is explicitly expressed as a mode superposition (See Appendix A.3, equation (A.215)). The strong inequality in (2.45) holds close to criticality. Suppose we were to restrict attention to the linear approximation, i.e., retain only terms that were linear in £. The linear approximation is e0+ai^h )2[ (K * ° -K *{K *0] = ° (2-46) Substituting in the full mode expansion for £(0) into this equation and equating coefficients of each harmonic on both sides of the resulting equation, leads to a set of equations, one for each mode amplitude £„. These equations are of the form: A„£n = 0, n = l , 2 , . . . (2.47) where the coefficient An coincides with an eigenvalue of the stability matrix, as shown in Appendix A.3. These equations have only trivial solutions, which indicates that the linear approxi-mation is not sufficiently informative. Thus, it is apparent that the few first nonlinear terms, which will lead to the interactions between harmonics are essential. It will be seen that quadratic and cubic terms are needed in the expansion and must be retained in the equation for £. As this equation is rather cumbersome, we leave its details to Appendix A.3 (equation (A.218)). The next step, as before, is to substitute the mode expansion of £(0) into the nonlinear equation and keep the first few nonlinear terms of the same degree of smallness. After some simplification, we equate coefficients of each harmonic on both sides of the equation Chapter 2. Models of self-alignment phenomena in cell biology 37 as before. This calculation leads to a system of equations (one equation for each mode). The system has a particular hierarchical structure which makes it straightforward to solve: |A2|6 + ^ 2 6 6 + ^ 6 6 = 0, |A3|6 + ^ 3 6 6 + ^ 6 6 = 0, (2-48) |A„|£n + -B„£l£n-1 + 5„6 (n+ i = 0. The linear growth rate of the first unstable mode Aj is the small parameter close to the bifurcation. Note that the first equation involves only £0, £1, £2, the second equation only £1, £2* £3, etc., so that they can be solved one at a time. The coefficients Bi, Bf, F, G are given explicitly in Appendix A.3, and are of order 1. Further, the fact that the total mass, M = C + P is constant and the assumption that |£„| <C |£„_i |, n = 2 , 3 , . . . lead to the estimate £0 ^ -At£l (2.49) where A\ is given in equation (A.219). The dominance of the leading mode over all other modes means that the terms whose coefficients are B\ in equations (2.48) above can be neglected compared to the other terms in each equation, further simplifying the task of solving the set of equations. We thus estimate the amplitudes of the stable modes as follows: in * f ^ i k - i - (2.50) We observe that this implies that the modes are strictly ranked in size, i.e., that |£n | ^C |£„_i | ,n = 2 , 3 , . . . . Plugging (2.49) and (2.50) for n=2 into the first equation of the system (2.48) we obtain a self-contained equation for £1 whose form is: A16 - D& = 0, (2.51) Chapter 2. Models of self-alignment phenomena in cell biology 38 (see Appendix A.3 for details). Careful estimates reveal that the coefficient D in this equation is positive unless the total mass M is very small ( M < C « ; but then C is very small and of the same order of magnitude as £, so that the method above then cannot be applied). Thus, the solution of (2.51) has the form: 6 = ± ^ (2-52) We now consider a and M as constants and introduce the idea of a small governing parameter, e, (which may be some combination of original parameters of the model). We assume that as this parameter crosses a critical value, ec, stability breaks. Very close to this critical value, i.e., when (e — ec) <C ec> Ai ~ (e — ec), and we can write (2.52) in the more revealing form 6 = ±kyj\t-€c\, (2.53) Further, using the estimate (2.50) gives us the additional result |£»| = fc»|e-ec|n/2, n = 2 , . . . (2.54) where k, kn are some coefficients of finite order. We have thus explicitly recast the amplitudes of all modes in terms of the governing parameter, e, and more specifically, in terms of its "distance" from the critical value at which bifurcation occurs. The bifurcation diagram based on these results appears in Figure 5a. It can be seen that the bifurcation is supercritical, implying a non-equilibrium phase transition of second order (see Friedrich and Haken (1989)). This means that as the parameter e increases past the critical value, the amplitude of the nonhomogeneity increases gradually from zero. This resembles the dynamics of growth of the order parameter in physical systems such as liquid crystals. So far we have assumed a single humped kernel. In the case of a double humped kernel, we have an expansion over the even harmonics only; the second harmonic is Chapter 2. Models of self-alignment phenomena in cell biology 39 then responsible for the stability break in general. In this case the reasoning and the mathematical derivations are entirely analogous to the case we have discussed in detail. The expansion is shown in detail in Appendix A.3 and the character of the bifurcations is the same. These results are in full agreement with the general scheme for symmetry break in systems with 0(2) symmetry (see Friedrich and Haken (1989); Busse (1987)). Biolog-ically, the results of the bifurcation analysis gives us the following pictures in the two dimensional case: In the case of the single-humped kernel, the steady state bifurcates into a distribution similar to cos(#), whereas in the double-humped case, it leads to a distribution similar to cos(29). This behaviour is illustrated in Figure 6a,b. 2.7 Bifurcation Analys is of Mode l I in 3 D The qualitative picture of the development of patterns in the 3D case following bifurcation is similar to the 2D case with one important exception. In the 3D case, the rotational symmetry group is 0(3) which leads to rotational and pattern degeneracy, i.e., several eigenfunctions correspond to the same eigenvalue. If the kernel is single-humped, then the first harmonic Y\ responsible for the stability break is doubly rotationally degenerate. There is no pattern degeneracy in this case. These statements are true if we ignore the two-fold degeneracy over the associated Legendre polynomials, P j ^ P j - 1 which differ in dependency over the angle 9, i.e., { cos# (2.55) sin0 where Y™ = P™(cos</>) cos(m#). These can be neglected since different linear superposi-tions of P * 1 just represent rotation about the axis of symmetry, thus leaving the essential character of the solution unchanged. Furthermore, this degeneracy does not represent Chapter 2. Models of self-alignment phenomena in cell biology 40 two different patterns. From group theoretic methods (Busse, 1987), it is known that the two harmonics Yj° and Y± represent essentially the same state, viewed from two different orientations. As before, we consider an expansion of C in terms of eigenfunctions, but now these are the Legendre polynomials, rather than simple trigonometric functions (See Appendix A.3). We are now concerned with the amplitudes of two unstable harmonics, namely P°(coscf)),Pl(coscj))cos(9) and we call these amplitudes Zi,z2, respectively. We further use the notation yn(t),n — 2 , . . . to denote amplitudes of the stable harmonics. In the steady state, Z\,z<i,yn are constants. We expect that the amplitudes of the stable modes again depend on the leading amplitudes, i.e. yn(t) = yn(zi(t),Z2(t))- Close to criticality (when the governing parameter approaches its critical value) the stable mode amplitudes will be strongly dominated by the leading mode, so that we have one of the inequalities, \yn\ < \zi\ or \yn\ < \z2\ or both: \yn\ < | ^ | i = 1,2, . . . To find out what actually takes place in a given situation, we must find the first terms in the expansion for yn(zi) and substitute them into the analogue of equation (2.43). (The form of this equation is the same but the leading partial derivative is replaced by the surface spherical Laplacian). This will lead to a pair of ordinary differential equations for the amplitudes z\, z2 from which we can find the stable stationary solutions. As in the 2D case, we keep terms up to third order, and arrive at an equation analogous to (A.218). Thus, the method is in all respects analogous to the one described in the previous section. Because of the greater complexity in the 3D case, the analysis is formidable, and we do not here include details. Fortunately, we can use the results of Busse, 1987 and references therein where the general system with 0(3) symmetry with quadratic nonlinearities was considered. Our equations are a particular case of this Chapter 2. Models of self-alignment phenomena in cell biology 41 general system. We describe these general results below. Let us chose (f> = 0 as the direction of the axis of symmetry of the pattern evolved. This means that the leading mode P®(cos<f>) = cos</> takes over. In this case, in the stable stationary state zi = 0 and we have the estimate 1/2 ~ z\- The relative sizes of stable modes are given by \yn\ <C I2/2I f ° r aU n > 3. The amplitude of the leading mode z-i obeys the equation \zi -dz^ = 0 (2.56) A is the growth rate of unstable modes in the linear approximation, d is a coefficient of finite order depending on C,a,M,Ki, AV Both A and d are positive. We observe the similarity of this result to that of equation (2.51). The important fact is that the damping nonlinearity which stops the growth of the unstable mode, and causes its saturation is a cubic one. The above equation has a stationary solution Zl = ± ( ^ / » . (2.57) The bifurcation is again supercritical, as shown in Figure 5a. From our choice of axis of symmetry, we have already noted that the leading mode Pj°(cos<^) = coscf> takes over. The pattern which starts to grow would therefore have the form shown in Fig 7a. This type of pattern has the following property: the number of cells whose orientation is <f> = | (i.e., the angular distribution at the equator of the spherical state-space) is unchanged, the number of cells at angles cf> < | is increased, and most greatly at the pole 0 = 0. Conversely, the density of cells at (j> > ~ diminishes, and particularly so at the pole 4> = IT. The case of the double humped kernel, is given in Appendix A.3. There is a difference Chapter 2. Models of self-alignment phenomena in cell biology 42 in that the leading mode possesses both orientational and pattern degeneracy. We con-sider the amplitudes xi, x2, X3 of the modes containing terms of the form P2°, P2, P2 in the mode expansion of C(<f>,6,t). It can be shown that as a result of intermode interaction, the pattern degeneracy is removed: all non-axisymmetric harmonics die out. Suppose we again chose <f> = 0 as the direction of the axis of symmetry of the pattern evolved. Then in the stable stationary state, the amplitudes x2 = £3 = 0. We can arrive at the following equation for the leading mode amplitude: \Xl + Plx2x - p2x31 = 0, (2.58) where pi and p2 — pi are positive constants given in equation (A.229) The stationary amplitude is then (A is small) -AM %i — < (2.59) P1/P2 The important feature of this solution is that, as in the case of a single-humped kernel, pattern degeneracy disappears. The harmonics Y2 and Y$ die out and the only leading harmonic with amplitude significantly larger than all other harmonics is P2°(cos <f>) = - ( 3 cos2 6 - 1). (2.60) Li The pattern evolved is axisymmetric. The expression (2.59) suggests that in this case, contrary to the 2D case (both single and double hump), and to the 3D single-humped case, we have a transcritical bifurcation as shown in Figure 5b. Examples in physical systems include liquid-solid or liquid-gas transitions. In physical language, this means that we have a nonequilibrium phase transition of the first rather than second order. This leads to a number of important qualitative conclusions. One of these conclusions is that even before the bifurcation (e < ec), there are values of the parameters (e'c < e < ec), such that the stable inhomogeneous pattern Chapter 2. Models of self-alignment phenomena in cell biology 43 can co-exist with the stable homogeneous distribution. Quantitative conclusions cannot be made from the diagram and formula (2.59) as the analysis is valid only at small amplitudes. This solution is characterized by the value x\ > 0, so that the concentration of the cells around the equator decreases from the value <f> « 55° to <f> « 125°, and the concentration of cells in the vicinity of the poles increases, see Figure 7b. Then the amplitude of this pattern is not small, and it does not grow smoothly, as predicted in the 2D case. Rather, the amplitude of the pattern jumps from zero to the values higher than xi quite suddenly. The process is analogous to nucleation, which is well-known in physics. That is, a small perturbation can lead to massive recruitment and thus grow in size abruptly, as parameters leading to bifurcation change a little. Our phenomenon, however, operates in angular rather than physical space. If t'c < e < ec, then, as a result of fluctuations which are larger as e approaches ec, too many cells gather around the poles, leaving the equator, and the system relaxes to the nonhomogeneous pattern. After the bifurcation another stable pattern evolves at which the amplitude x\ « — X/pi is small and negative, so that the concentration of the cells decreases around the poles and increases around the equator. The results of our nonlinear analysis are in qualitative agreement with the general mathematical results of pattern formation in systems with 0(2) and 0(3) symmetries (Friedrich and Haken (1987); Busse (1987)) and well-known physical results on the char-acter of phase transitions in ferro-electrics and liquid crystals (de Gennes (1974); Chan-drasekhar (1977)). These results of nonlinear analysis apply to Models II and III. Chapter 2. Models of self-alignment phenomena in cell biology 44 2.8 Discuss ion According to our linear stability analysis and bifurcation analysis, all three models in this chapter lead to qualitatively similar results. This signifies some robustness in the modelling approaches. Not only are the results independent of the detailed assumptions about turning rate kernels, but also the types of forces or effects leading to turning do not significantly influence the behaviour. In short, the result of this chapter is that as interactions effectively increase (with the growth of the density, for example), the homogeneous angular distribution becomes unstable, and mild inhomogeneities evolve. These inhomogeneities signify the onset of the growing angular order. The type of the inhomogeneities is defined by the symmetry properties of the kernels. We would like to remark on the similarity of cell alignment to the alignment phe-nomena that is described in the physics literature: that of liquid crystals (de Gennes, 1974). The particular case of nematic liquid crystals, in which the centers of mass of the molecules have no particular order, but in which there is some order in the directions of orientation of molecules is an apt analogy to the dense cultures of fibroblasts that ex-hibit patch alignment (Elsdale, 1972). Molecules in such structures can undergo random motion and tumbling, and they interact by electrostatic attraction or repulsion. In cultures of fibroblasts, the random turning rates of the cells are analogous to the tumbling of molecules, resulting in a kind of angular diffusion. The cells however, are living units, with essentially "infinite resources of energy" on which to draw. Their inter-actions cannot be easily summarized with simple physics. The alignment of populations of cells is not an outcome of the shapes of the cells, but of the complex membrane and cellular cytoskeleton, and the dynamic response to contact with another cell. In work dating back to Onsager (1949) and Zwanzig (1963), the case of hard rod-like molecules which do not overlap was studied using thermodynamic principles. It was Chapter 2. Models of self-alignment phenomena in cell biology 45 shown that thermodynamic considerations of entropy alone, without forces of attraction between molecules could account for the long-range orientational order in these liquid crystals. Several rigorous models in the physics literature predict orientational phase transitions of the second order for nematics, and these are analogous to the ones we observed in our models. A key difference in the methods of approach and the tools used for analysis of these physical phenomena must be emphasized. Models in physics are traditionally based on minimization of a free energy functional, which mandates that the system studied is being investigated close to thermodynamic equilibrium. However, in our approach the interest is on the dynamic process itself, far away from such equilibrium. The reason that we abandon the traditional physics approach when dealing with these biological systems is that a meaningful definition of free energy cannot be derived from first principles in highly non-equilibrium open systems such as living cells. It would be possible to formally define a free-energy functional, perhaps, as has been done for some open non-equilibrium systems, but this approach may be artificial. This alternate energy approach is explored in Murray (1989), Section 9.6. Chapter 3 Limit ing case of nearly complete al ignment 3.1 Introduct ion: the Peak Ansatz The results of this chapter are original, except for formulae in the introduction and results of Section 3.6. The latter were obtained by G. B. Ermentrout and appeared in a joint paper (Mogilner et al., 1995b). I describe them here for the completeness of a general picture. Also, numerical experiments described in Section 3.7 were carried out jointly using the program developed by L. Edelstein-Keshet. In this chapter I am particularly interested in sharp peak-like solutions that develop when the effects of dispersal are very small. Such solutions are known to occur in a number of systems of physical and biological relevance. Examples include: peaks of cell density in chemotaxis equations (Grinrod, et al., 1989), in cell-contact models (Edelstein-Keshet and Ermentrout, 1990), and in models for individual aggregation in population dynamics (Grinrod, 1988). Other applications include the analysis of models for the distribution of traits (for example, dominance) in animals (Jaeger and Segel, 1992). First, the technique, which we call the Peak Ansatz , will be illustrated on an explic-itly solvable system from the literature, and then this method will be applied to studying behavior of the models for angular distributions described in the previous chapter. The following example comes from a model for the dispersal of a population of insects with density distribution u(x,t) given in (Grinrod, 1991): ut = euxx - (uw)x (3.61) 46 Chapter 3. Limiting case of nearly complete alignment 47 where foo rx w = u(y)dy - / u(y)dy (3.62) Jx J—oo Here w is the swarming velocity, and is positive if the total number ahead of an individual is higher than the number behind it. (If the number ahead is lower, the velocity is negative). The parameter e is a measure of the random dispersal of the swarm. The model is somewhat unrealistic as it assumes an infinite viewing horizon to the front and back of an individual. However, it is attractive mathematically, as it is an explicitly solvable case. The steady state solution (given in Grinrod (1991)) is: M2 M u(x) — —— sech (~—(x — XQ)) (3.63) 8e 4e where /oo u(y)dy (3.64) -oo is the total number of insects in the population. For us the most interesting case is the case of e approaching zero, when the solution (3.63) represents a <£-like peak with width of order e. In the limit e = 0, the shape of the peak is formally a 5 function itself. The location of the center of the peak, XQ is arbitrary due to the fact that space is homogeneous. In this particular model, the unrealistic but simple nature of w permits the explicit solution given above to be computed. However, in many more realistic models, no such explicit solution can be found, and then the peak ansatz can reveal the behavior of solutions for small dispersal rates. Our approach to this problem will be as follows: We expect to find the solution of the truncated equation (e = 0) ut = -(uw)x, (3.65) in the form of a superposition of a few delta-like peaks, the heights of which are constant. The locations of the peaks change due to interaction between individuals. The only Chapter 3. Limiting case of nearly complete alignment 48 stationary stable configuration would be that of one peak. If one now wants to take diffusion into account, and this diffusion is slow, one can consider it as a small stochastic perturbation of the quasi-deterministic processes described above. Mathematically, one would have a singular perturbation problem. (The character of the problem changes if the small parameter e is set equal to zero.) This situation is similar to the Fokker-Planck equation with small diffusion term described in Gardiner (1985). As in his treatment, we introduce a scaled variable and try a perturbation expansion of the solution. One of the approaches to this type of problem is to assume solutions of the form u(x) — exp(—</>(x)/e), (3.66) where e is a small parameter, and then look for an expansion of the form # r ) = S » » ( x ) . (3.67) (See Gardiner (1985)). Our approach is similar to this, but rather than dwell on the mathematical technicalities, I essentially consider only the first approximation in this expansion. The idea illustrated here will be applied to studying the behaviour of three angular-distribution models which were derived and introduced in the previous chapter. In three models, introduced in Chapter 2, there are two types of elementary processes leading to turning of the individuals: The first is random rotational diffusion, while the second is deterministic turning caused by direct interaction. A special case occurs when rotational diffusion is very slow compared to the interaction-mediated turning. In this case, the random element of the models can be essentially neglected. Models I and III have the following important feature in common: When rotational diffusion is omitted, individuals can only "be attracted" to some angle by sticking to other individuals or clusters already oriented at that angle. Thus, if there are no individuals at some angle, Chapter 3. Limiting case of nearly complete alignment 49 none will appear there. This suggests applying the peak ansatz and reducing the system of integro-partial differential equations of the models to ordinary differential equations. In the limiting case of e = 0 in Model II, if a number of cells are concentrated at the same angle, the angular velocity of each individual in this group is completely prescribed by deterministic terms in the model. Moreover, individuals at a given angle move collectively with the same velocity. This can be represented as convection of a number of < -^like peaks in the angular distribution. Contrary to the case of Models I and III, in Model II the location of the peaks changes but the heights remain constant. Each peak has a deterministic trajectory prescribed by a kinematic equation (which is Liouville-like). There are a number of stable static equilibria for the ensemble of <£-like peaks. In a way, this case reduces to a problem in nonlinear classical mechanics. Again this limiting case allows analytic treatment. 3.2 Appl icat ion of the Peak Ansatz to Mode l I Let us start with equations (2.20) of Model I and explore the behaviour of the solutions to these equations in the limit as the parameter e gets small. This means that the angular diffusion of the population becomes insignificant compared to the rates of the other processes. Note that the absence of diffusion does not mean that there is no transfer of mass from one angle to another, since the convolution terms still represent reorientation of cells due to interactions. (For small rotational diffusion cells tend to maintain their directions of motion until they make contact with another cell.) This suggests a solution in the form of a number of $-like peaks in the angular distri-bution, the location of which is constant and the heights of which change due to exchange of individuals. The heights of the peaks are governed by a system of nonlinear ordinary differential equations (ODEs). Chapter 3. Limiting case of nearly complete alignment 50 The stability properties of these peaks can be studied. It is convenient to rewrite the model in terms of the total number at a given angle 0, which is M(6, t) = C(0, t) + P(0, t). (3.68) Plugging the expression P = (M — C) into (2.20) and assuming e = 0, the model can be reformulated in a more convenient form: °-g = M{K*C)-C(K*M)i g = - C ( a + ( A ' * M ) ) + aM. (Note that while not evident from this form of the model, the positivity of M and C is implied by the original model.) The second equation can be interpreted in the following way: the rate of change of the total mass of cells at a given orientation is the sum of effects of free cells attracted to that direction minus the mass of free cells that are attracted to other directions. The expression K * M represents the influence of the mass distr ibution on the direction 6. Since it appears repeatedly in these equations and their analysis, it is convenient to use the notation F(O) = (K*M)(0). (3.70) Setting the time derivatives to zero, we arrive at the steady state equations: ~C{a + (K * M)) + aM = 0, M(K*C) = C(K*M). From equations (3.71), we can observe that (3.71) C = aM = aM (3 72) (a + (K*M)) (a + F)' { ' ' Further, inserting this result into the equation for M, we get ""•<^T>-"&g- <-» Chapter 3. Limiting case of nearly complete alignment 51 The only possible solutions to this equation are those satisfying pointwise M = 0, (3.74) or K-, M _&*M) Here only the class of peak-like solutions will be explored: ( o , o#{0i,...,en}, M = ^ (3.76) Note that I am not claiming that this solution is unique, as there may also be solutions with a finite support. However, in the case of this model of a finite collection of objects, such solutions are of lesser interest. As will be shown by numerical simulations, these peak-like solutions evolve from a variety of initial data. Recall that this type of solution is here refered to as the Peak Ansatz . That is, the mass is concentrated in a series of k peaks (8 functions) at the angles {0\,..., 0n}. Biologically, this solution can be interpreted as the existence of a discrete set of direc-tions {9i, 02,... 0n}, n — \,2,... along which cells are oriented. The masses concentrated at these angles, namely {Mi, M2, • • • Mn} satisfy a set of n equations: (K * M){0i) = constant, 0{ e {0U 02,... 0n}. (3.77) 3.2.1 Reduct ion of Mode l I t o discrete equations governing the peaks In this section the Peak Ansatz will be applied to reducing Model I from a set of integro-partial differential equations to a set of ordinary differential equations for the masses concentrated at a number of peaks. I now will consider individually the cases of one, two, and more than two 8-like peaks. An eventual goal is to investigate the stability of such peaks. Note that understanding stability to an arbitrary perturbation is a challenging Chapter 3. Limiting case of nearly complete alignment 52 problem not at tempted here. Rather, I formulate and study the equations governing interactions between these discrete peaks, and determine whether one or more peaks can persist. I develop here a short-hand notation for describing various special cases. N o t a t i o n for the peak-like solutions In the equations below the following short-hand notation is used for the probability that a cell at direction 9{ turns to direction Of Kij = K{pi - Oj). (3.78) Observe that since K{0\ — 02) = K{02 — #i) it follows that K,j = Kf. The convolution in equation (3.77) is now a finite sum, rather than an integral, but the interpretation is the same: the influence of the distribution is the same for each value of #;. Thus, the steady state satisfies £ KuMj = J2 K2iMj = ... = £ Kn]My (3.79) j=\ j=i j = i The conservation of total mass must also be maintained, so that £ Mj = M. (3.80) These steady state equations have a unique solution in the set of real numbers, but clearly, for biological applications, only nonnegative values of the M / s are acceptable. This system does not have a unique solution if the matrix {G!j};j=i...n is singular, where Gij = K^ - Knj, i = l...,ra - 1, j - l...n; Gnj = 1, j = l...n (3.81) The elements of this matrix are defined on an (n-1) dimensional parameter space of (n-1) angles between the peaks. This matrix is singular on some manifold of co-dimension 1 in the parameter space (Arnold, 1978), so, in a generic case, steady state equations (3.79-80) have a unique solution. Chapter 3. Limiting case of nearly complete alignment 53 The amplitudes of the peaks, {Mi, M2,... Mn} are fully determined by the set of directions {#1, #2, . • .9n} given a particular choice of the kernel function K{9), as the values of the coefficients Kj are defined by this set of 9 values. This means that not all sets {#i, 92,... 9n} have meaningful solutions. It can happen that an arbitrary set of directions leads to a system of algebraic equations having one or more negative values of Mj in its solution. This set must then be rejected as unbiological. Given a nonnegative solution, {Mi, M2,. •. Mn}, the masses of the free cells, {Ci, C2,... are given by the expression C3 = -^-M3, F0 = J£KijMi. (3.82) a + -To j=i In order to understand stability of the steady state in the general case, in the sections below we will consider sequentially the behaviour of one, two, and three peaks. I will restrict attention to a limited class of perturbations consisting of a finite number of peaks. The analysis of the general case is considerably more challenging. 3.2.2 Stabil i ty for the case of two compet ing peaks In this case all material is concentrated at two peaks located at angles 9i and 02, and the amplitudes (Mi,Ci) and (M2 ,C2) define the distributions. According to the definition of K^, cells having the same orientation will be represented by interactions Kn = K22 = K0, (3.83) and the interaction of cells in distinct clusters is given by: K12 = K21 = Ki. (3.84) Recall the assumption that A'(0) represents interactions that are stronger than at other relative orientations. This means that K0 > K^ i ± j . (3.85) Chapter 3. Limiting case of nearly complete alignment 54 Then by equation (3.79), K0M! + KXM2 = K1Ml + K0M2. (3.86) The equation above leads to *, = #, = £. ft = ft = - j ^ , (3.87) where K = (KQ + K\) and M is the total mass. Below we use the symbols C and M to denote steady state solutions. Now time-dependent peaks at fixed angles are considered. I start with the properties of the stationary solution given in (3.87), namely the case of two peaks. Since mass is concentrated only at two discrete angles, we can consider the ordinary differential equations that describe these masses. It can be easily seen that M = (Mi + M2) = 0 (where the traditional dot notation is used here for time derivatives.) This can be put into a more convenient form by defining the difference of masses m, as follows: Ml-M2 m = . (3.88) The dynamics of the variables, m, Ci, C2 is given by the system m =K1[f(C2-C1) + m(C1 + C2)], Cx =2M + am-(a + £f)C1 + (K1-Ko)mCu (3-89) C2 = ^ + am-(a + ^M)C2 + (K0 - Ky)mC2, I look for the solutions of this system in the form C\ = C1 + C1, C2 = C2 + c2, m, where m,Ci,C2 are small perturbations from the steady state. In Appendix B.l the linearized system of equations for the variables m,ci,c2 will be derived, and linear stability will be investigated. It is shown there that small deviations from the stationary steady state grow exponentially and that therefore two peaks are unstable. Chapter 3. Limiting case of nearly complete alignment 55 3.2.3 Stabil i ty for the case of a single peak Here I formally consider two directions in order to investigate the stability of one peak. In the stationary state, I consider all cells concentrated along the first direction and none of them along the second direction. As a result of a small perturbation, a small peak may appear in the second direction. When the mass is entirely concentrated in one peak (all cells are lined up along a single direction), then M1 = M, M2 = 0, &= "^ .., C2 = 0. (3.90) a + KQM The stability of one peak is to be considered. I look for solutions of the system (3.89) in the form C\ = Ci + ci, C2 = C2+c2, m = fh+m, where m, Ci, c2 are small perturbations of the steady state (3.90) (in this case the steady state value of m = -^-)- In Appendix B.2 the linearized system of equations for the variables m, ci, c2 is considered in order to investigate linear stability. The conclusion reached there is that one peak is stable to a perturbation consisting of one other small competing peak. 3.2 .4 T h e case of three interacting peaks Two possible cases are considered, one in which the three peaks are equally separated from each other by angles of 27r/3, and a second case in which there are three uneven peaks. (a) Equal ly spaced peaks In this case, 62 — #1 = #3 — 02 = 9\ — 03 = 27r/3. We will denote the values of the kernel as A'(0) = K0, K(±2TT/3) = Kx. Then it follows from (3.79-80) that a steady solution exists provided A'0M1 + A'1(M2 + M3) = K0M2 + K1(M1 + M3) = K0M3 + K1{M1 + M2). Chapter 3. Limiting case of nearly complete alignment 56 As KQ > Ki, the stationary solution is Mi = M2 = M3 = M / 3 . (3.92) I am not interested here in solutions having some of the Mj-'s zero, since these cases have already been covered under the previous analysis of two and one peaks. Then the masses at the three peaks satisfy a system of six equations which are shown in Appendix B.3. The stability calculation there demonstrates that three equally spaced peaks are unstable. As expected, the analysis of this problem is more taxing than those of the previous cases of one and two peaks. (b) Unequal ly spaced peaks In this case #2 — 9\ ^ #3 — 92 7^ #1 — 93 7^ 92 — #i- Then there are three peaks located at some angles which are not equally spaced. The sizes of the peaks may also be different. This case in which the peaks are unequal and/or are located at three arbitrary angles is analytically forbidding, and will be explored with a limited numerical experiment. The problem of investigating the stability of three or more peaks in this model appears to be extremely difficult and may or may not be analytically tractable. This problem has a full analytic solution in the case of Models II and III. But from the previous discussion it would appear that the only stable situation is a single peak, i.e., a stationary solution in which all cells are lined up along a single direction. The equations of this system are analogous to those of case (a), but allowing for differences in the values of K(O),K(0i — 92),K(9i — 93),K(92 — 93), since the angles need not be equally spaced. As a specific example, the software program Mathematica for certain parameter values was used to investigate the stability matrix. The exact values and the resulting stability matrix and eigenvalues in that example case are shown in Appendix B.3. Two of these eigenvalues are positive, so that this situation of three uneven peaks is unstable. One of the eigenvalues is zero since M = Mi + M2 + M3 = const. Chapter 3. Limiting case of nearly complete alignment 57 3.3 T h e form of the peaks in Mode l I (weak angular diffusion) To find out more about the form of peaks in Model I before they become infinitely sharp (when e attains the limit e = 0) consider the influence of a small but non-zero value for the random turning rate, e ^ 0 in equations (2.43-44) for the distribution of free cells. The presence of a small parameter as a coefficient of the second derivative term leads to a singular perturbation problem. (See, for example, Chang and Howes (1984); Smith (1987); O'Malley (1974)). Clearly, the presence of the Laplacian operator will cause the sharpness of the peaks to be smoothed somewhat, depending on the value of e. The distribution of bound cells can also be easily found with the help of expression (2.42). (In particular, it is found that the width of the peak of P is of the same order of magnitude as the one with the free cells.) In the limiting case e = 0, it was found that a stable steady state exists corresponding to all cells concentrated in one J-like peak with zero width. In the presence of small diffusion, this problem can be analyzed using singular perturbation theory, due to the presence of the small parameter e in the highest derivative term in equation (2.43-44). The fact that this is a nonlinear integro-differential equation makes this a difficult problem. The unperturbed solution is taken to be a <£-like peak which is non-vanishing in the limit e —> 0. The solutions are asymptotically zero except at regularly spaced points (one point in our case). In the neighborhood of such a point the solution has a spike. Without loss of generality, assume a peak at 9 = 0 and (7 = 1. It will be seen that a shape of the peak can be approximated by the form A(e)exp(—02/a(e)). Then \[G represents the width of the peak. I am interested in how a depends on e. Since C(9) is asymptotically zero away from 9 = 0, the behaviour of the function V(C) is of interest in the vicinity of 9 = 0. Chapter 3. Limiting case of nearly complete alignment 58 Rescale the angular variable x = -jy^ -, and look for the solution of equation (2.43) in the form S(x) = e1/,4C where S obeys the following equation: e 1 / 2 | ^ + V(e^x)S = 0. (3.93) (With a slight abuse of notations we use the form V{tl'Ax) for the function V{C{tl'4x)).) Assume that S can be expanded in the following asymptotic series: S(x) = ^=0en'ASn(x). (3.94) and that V(e1/4x) = Z™=2en/4Vn(x). (3.95) I find only the first terms in the asymptotic series for S and V. The equation for the first approximation has the form fgi + V2(x)S0 = 0, V2(x) = ax- a2x\ e^x « 1, ai = bA2, a2 = | , 6 = - a ^ p > 0 , A2 = fZox2S0(x)dx. (The form of the function V2(x) was obtained from plugging the expansion (3.94) into the expression (2.43). Then the kernel K(e1'4(x — x')) was expanded in Taylor series. After cumbersome calculations two terms of order e1/2 were obtained). Equation (3.96) is the well-known equation for the harmonic oscillator in quantum mechanics. (See Morse and Feshbach, 1953). There is a unique strictly positive solution to this equation which decreases at infinity. This solution occurs under the condition ct\ = yjci^,. This condition is equivalent to the expression A2 — l/\/2b. Under this condition, the normalized solution is s° = ^ ) 1 / 4 e ^ - y f x 2 ) - (3-97) Chapter 3. Limiting case of nearly complete alignment 59 Calculating A2, using the definition (3.96) and the form of the solution (3.97), we find that . 1 / 2 1 A2 = 2\lb iM' Thus the condition a j = ^Jcl^, is fulfilled, and so the solution fulfils a consistency condi-tion. The corresponding form of the peak in the first approximation has the form: C° = ^( A ) 1 / 4 e ^ ( - ^ 2 ) - (3-98) This means that the width of the peak in the limit of slow diffusion in Model I is small, of order e1/4. 3.4 Appl icat ion of the Peak Ansatz to Mode l II N o t e : In this section for the sake of convenience, we write the M o d e l II in the form: F)C 82C f) -ai = 'W-wMw,c))< (3'99) where w * c = I w(e - e')C($', t)de'. (3.100) We assume that W{6) is an odd function, and that W(0) = W(ir) = 0, W{6) < 0 for 0 < 0 < 7T, and W(6) > 0 if -7T < 0 < 0. In this section, we consider Model II in the limiting case of zero rotational diffusion. The equation for the dynamics of the continuous angular distribution has the form ^ - - ^ 1 (3 101) dt ~ 80 ' ( j where v = W*C. (3.102) Chapter 3. Limiting case of nearly complete alignment 60 From these equations, it can be seen that the angular dynamics become purely deter-ministic in this case. This fact suggests that the motion of each individual object can be followed. This is the so-called Lagrangian approach, in the traditional terminology of theoretical mechanics in contrast to the Eulerian approach of following behavior at a given point used in the rest of this Thesis. Mathematically, the nonlinear equation (3.101) is hyperbolic, and the trajectories of individual objects are given by the equa-tions of the characteristic curves. If the equations of these characteristics are solved, a link can be made between the two approaches. Unfortunately, in general, nonlinear equa-tions for characteristics can not be solved in closed form (Logan, 1994). For this reason, we restrict attention to a special case whose qualitative nature can be fully described. Let the vector describing the state of the system at a given time, t, be n(t) = {0u...en}. Here, 9i(t) is the orientation of the ith object at time t. The number n of objects in the system is conserved in time (n ^> 1). The dynamics of the ith object is given by the definition of its angular velocity: ^ = "MO)- (3-103) To ensure that the former Eulerian definition of density reduces to the Lagrangian defi-nition, it is necessary to define the angular velocity as follows: i,(0) = r W{0 - 0')C(6')d6', J — -K where C(9) is formally defined as: C{6) = ^=l8{6-ei). Then the angular velocity can be written in the form: v(e) = X]=1W(9-0j). (3.104) Chapter 3. Limiting case of nearly complete alignment 61 We introduce the "potential function" V(9), V(6) = - Jw{6)d0. (3.105) (The constant of integration is irrelevant, as potentials are always defined up to some additive constant.) We also define the Lyapunov function V(tt): v(n) = Eiij)v(ei-0j). (3.106) (Here the summation is taken over all pairs of objects.) Now the dynamics of the objects is governed by the following equations: d9j _ dV(tt) dt d9i " These equations can be written in the vector form: The advantage of this approach is clear from the form of the last equation. This is now a gradient system (Gross and Hohenberg, 1993). The behaviour of such systems is much simpler than that of general non-linear dynamics. If an attractor ftQ can be found such that the Lyapunov function V has a global minimum at fl = Qo, then the dynamics of the system is just relaxation towards this attractor. With the definition of the interaction kernel W{9) at the beginning of this section, the corresponding potential function V(9), (3.105) has a global minimum at 9 = 0. Note the analogy with the following physical description: the interaction kernel W corresponds to the forces between two objects. According to the definition of the kernel, these forces are attractive. The potential V thus has a global minimum when the distance between two objects is zero. The Lyapunov function V is the total potential energy of the system. Chapter 3. Limiting case of nearly complete alignment 62 Obviously, in the present case this total energy has a global minimum when all objects have the same orientation. The global attractor for our definition of the Lyapunov function fl0 = {9,9,...0} where 0 is an orientation of complete alignment (arbitrary due to rotational invariance of the system). The above argument supports the assertion that the objects converge to a unique globally stable state of complete alignment. Note that another locally stable stationary angular configuration may exist. The biological meaning of this fact is that in nature there may be quasi-stationary angular patterns different from total alignment. 3.5 T h e shape of peaks in Mode l II (weak angular diffusion) Consider now the situation when one globally stable peak is present and angular diffusion is weak, but not absent. This situation again poses a singular perturbation problem. As in the case of the first model, I look for a solution in the form of a single sharp peak at 0 = 0. Integrating equation (3.99) once leads to the first-order ODE: 3C s— = vC, v = W*C. (3.108) Ov (The constant of integration has the meaning of flux, and it is easy to show that this con-stant is equal to zero for periodic solutions on the circle due to the rotational invariance). Rescale the angular variable: x = -4- and look for a solution in the form S = y/eC, where S obeys the equation dS _ v_Uexls = Q Ox y/e The region of interest for the solution of this equation is y/ex << 1 , 6 << 1. We assume Chapter 3. Limiting case of nearly complete alignment 63 that the functions v and S have the following asymptotic expansions: v(V~ex) = - £ ~ l U n e ? x n , S = £ ~ 0 e ? S n . Here I find the first approximation 5*0 which obeys the equation S'0 = -v&So, (3.110) with the solution 2 So ~ exp(^j^). (3.111) The first term in the expansion for the function v can be found as follows: v-^± = I r W(y^(x - x'))S(y/£x')dx' y/S £ J — oo "I /-oo = - / W(y/e(x - x'))S0(x')dx' + 0{y/l) £ J—oo 1 r°° = ~^W'(0) / (X - x')SQ{x')dx' + 0{y/t) = W'(0)Cx. y/£ J—oo Here I used the facts that W(0) = 0, that / xSo(x)dx = 0 and the definition of C: C = fC(x)dx. Thus, Vl = -W'(0)C > 0. (3.112) (Recall that W'(0) < 0). The peak solution has the form: V 7re 2e We have shown that in the limit of slow diffusion, the width of the peak is small, of order e1/2. Chapter 3. Limiting case of nearly complete alignment 64 3.6 Appl icat ion of the Peak Ansatz to Mode l III I now briefly discuss Model III when angular diffusion is absent (e = 0). The results of this section were obtained by G. B. Ermentrout and appeared in our joint paper (Mogilner et al., 1995b). If there are initially n infinitely narrow peaks in the angular distribution, then the redistribution of individuals over various angles (i.e., the transitions of the individuals from one peak to another) is due to the integral term in the equation r)C — = C(Q(C)*C). (3.114) This equation describes a set of discrete peaks (say n peaks) at constant angles, {9\,..., 9n} with {ci(t),..., cn(t)} individuals in each peak. None can appear in other directions due to the form of the equations. This model can be treated analytically by considering a discrete set of equations for the sizes of the peaks. Then the continuous convolution in (3.114) will be represented by a discrete sum, leading to an equation of the form dci n -E7 = ciJ2 L(°i ~ ci)GHcJ' (3.115) m j=i Equations (3.115) (one equation for each peak) state that the growth and decay rates of a peak c2 are proportional to its size. Whether this peak actually grows or decays depends on competition with all other peaks. This model was proven to have the following properties: 1. Any steady state of (3.115) is an equal subdivision of objects among n peaks. This means that ct- = 0 or c,- = M/n. 2. The steady state in which the population is concentrated in one peak is stable. 3. All other steady states are unstable. Chapter 3. Limiting case of nearly complete alignment 65 Thus, in this third model, in the case of no angular diffusion, the only distribution which is stable is one in which all individuals are in a single cluster, having the same orientation. 3.7 Numerica l exper iments for Model I In this section I compare the predictions of the analysis with a few numerical experiments, in which equations (2.1) were simulated in the 2D case, starting from an initial situation in which a number of peaks were present. (A set of three roughly equally spaced peaks, shown dashed in Figure 8, was used). Simulations were carried out by methods described in EKE (1990). The corresponding program was developed by L. Edelstein-Keshet. The results give a good qualitative description of the model. I started with a kernel which had a single hump and was strictly positive for all 0, K{6) = 1.5 + cos 6. This means that there are interactions between cells at all angles. Figure 8 shows how at weak diffusion, a single narrow peak evolves from a set of three equally spaced, small peaks. (Resultant peak can be anywhere depending on initial conditions. If the initial conditions are homogeneous, the "computer noise" defines the location of the peak.) Next, in Figure 9, the results of a simulation using a kernel with one hump which is nonzero for a finite range, smaller than 2TT are shown. The form of the kernel is K(0) c o s ^ , if - 1 . 7 4 < 0 < 1.74, un (3.116) 0, otherwise. (The range of influence was taken to be 100°.) It can be seen that two out of the initial three small peaks are absorbed into one large peak, and the remaining small peak separately grows into a smaller final peak. In this degenerate situation though two Chapter 3. Limiting case of nearly complete alignment 66 unequal peaks are not stationary, their rate of equalization becomes small. Therefore, in practice, two peaks of unequal size can "coexist". In Figure 10, the range of the kernel was decreased to 55°.The form of the kernel in this case was: K{9) = < cos-|r , if - 0 . 9 6 <9 < 0.96, •61 (3.117) 0, otherwise. In this case, three meta-stable peaks grow in positions identical to those of the three small peaks. The explanation for this phenomenon is that due to the narrow kernel, the peaks have a minimal interaction with one another. This was not true in the previous simulation, where the tails of the peaks were separated by an angle smaller than the effective range of interaction. Note that the peaks are not of the same height. In Figure 11, I explore the same situation as in Figure 10, but with a larger diffusion coefficient (e = 1). The initial peaks apparently start to communicate as a result of increased diffusion. Finally, a single peak evolves. In Figure 12 a double-humped kernel K(0) = 1.5 + cos(2#) is used. In the case of small diffusion two unequal metastable peaks evolve. Though these peaks are non-stationary, the diffusion is so small that exponentially small growth rate does provide a sufficiently strong equalizing effect on the two peaks, in real time. In Figure 13 I have the same case as in Figure 12, but with much larger diffusion. Evolution of two equal peaks at mutual angle IT can be clearly seen, as will be mentioned in the discussion. Note that diffusion in this case is much larger than in the previous cases, providing effective communication between the peaks. It is seen that the evolved peaks are slightly uneven as the relaxation time is very long. The results indicate that if the kernel is effective only at a finite range (as in the case of Figures 9 and 10) then in Chapter 3. Limiting case of nearly complete alignment 67 reality, a few metastable peaks will evolve at small diffusion. The lifetime of the peaks will be larger than the timescale of most realistic experiments. In real situations, this implies the possible existence of a few different directions of alignment. 3.8 Discuss ion The cases considered in this chapter all pertain to phenomena in which individuals achieve total alignment with one another. This stems from the fact that interaction terms domi-nate over random rotation. The alignment is an active one, stemming from forces, rather than a result of packing constraints due to crowding. This does not imply that all in-dividuals align in the same direction. They can be distributed among a few directions, depending on symmetry features of the interaction terms. The actual directions along which alignment occurs are arbitrary (unless we incorporate some bias due to external influences), but the relative angles of separation of two or more peaks result from the interactions. I explored the three distinct models for alignment in which turning was instantaneous (Model I, III) and gradual (Model II). In Models I and III a number of peaks of constant location and variable heights was found but the only stable configuration consists of a single peak. This fact is established definitively for Model III and in a few simple cases in Model I. For Model II, peaks of constant heights appear, and gradually drift towards each other, and finally converge into a single peak. To place these results into a biological context, recall that the angular densities represent concentrations of filaments (e.g., actin) or cells (e.g., fibroblasts) at various orientations. In Models I and III, these objects cluster at several directions, but small fragments are exchanged, so that one or another of the clusters grow and eventually win over the others. In Model II, clusters are gradually "pulled" by one another, so that , eventually, all have the same orientation. Chapter 3. Limiting case of nearly complete alignment 68 It follows from the qualitative theoretical reasoning in this chapter that the case of weak angular diffusion is characterized by a hierarchy of times in the dynamics of alignment. One may speculate that first, on a relatively short timescale, an "initial" configuration of peaks evolves from the initial distribution. Then these peaks undergo slow evolution, changing their location and heights on an exponentially large timescale (i.e., exp( l /e)) . The proof of this assertion is one of the goals of future research. The important fact for applications is that this large timescale can be comparable, or even larger, than the characteristic biological times. In this case, the system just does not have "enough time" to reach globally stable configuration. Thus, it may happen that the system aligns, but that the type of alignment seems to depend on initial conditions. There is an analogue to this scenario in small noise perturbation theories (Gardiner, 1985). There, also, a hierarchy of relaxation times in the system appears. The connection is that the Fokker-Planck equation with small diffusion is considered there, and the smallness of the diffusion results in this hierarchy. The nonlinearity of the problem results in great difficulty in performing a full treat-ment of the singular perturbation problem. Even the enumeration of all possible equilibria in the non-perturbed model is challenging. For this reason, I do not try to construct the full perturbational expansion and consider only the first approximation. I am not aware of theorems on the existence and uniqueness of solutions (or of Lyapunov functions) for Models I and III. This nonlinearity makes the problem of the stability of the peak ansatz solution prob-lematic. It is not clear how to prove global stability. Moreover, analytically only the local stability towards the shifts of the peaks, and their splitting can be checked. Theoretically, other perturbational modes may exist, breaking the peak's stability. However, there is a Chapter 3. Limiting case of nearly complete alignment 69 posteriori strong evidence for the stability of the peak solutions from numerical exper-iments. Further, there is no proof that the peak solutions are in a basin of attraction (with respect to most initial conditions), but again, numerical experiments seem to point to this fact. Now I briefly describe the expected results on first, parallel alignment due to a double-humped kernel and, second, orthogonal alignment that occurs in Model I for kernels representative of the alignment of actin filaments in the cytoskeleton of the cell (see Civelekoglu and Edelstein-Keshet, 1994). I make the following conjectures: 1. In the case of parallel alignment, the only stable situation is two equal peaks at relative angle of n in 2D and similarly antiparallel in 3D. 2. In the case of the orthogonal alignment I expect that an even, nonnegative dou-ble humped kernel having humps at —7r/2 and 7r/2 leads to orthogonal alignment. (Alignment occurs along angles where the kernel has maxima.) Far from criti-cality, where e —> 0 (diffusion is weak), four smooth bumps in the distribution of cells generate four equal and mutually orthogonal narrow peaks, located at 0 = {O,±ir/2t TT}. 3. In the 3D case of the orthogonal alignment, in the limiting case of weak diffusion, e —> 0, six equal, mutually orthogonal narrow peaks evolve from the corresponding bumps. In the most convenient case, when one of the peaks is located at the pole, their location is given by the set of coordinates {4> = 0}, {</> = 7r}, {</> = n/2; 6 = 0 , T T , - T T / 2 , TT/2}-4. It is plausible that the stable configurations of peaks in 2D and 3D situations described above in the weak diffusion limit are the only stable ones. Chapter 3. Limiting case of nearly complete alignment 70 A number of questions remain open for further investigation. One is the problem of generalizing the results to the 3D case. Perhaps, as in the 2D case, a single peak will be stable and its width will be given by the 2D expressions. In the cases of the double humped and the orthogonal kernel, other than proving the results rigorously, the problem of the shape of the evolved peaks is open. The global stability of a single peak in Model II for zero diffusion has been proved. It would be interesting to find out whether this stability is present in the case of small diffusion. Finally, more realistic spatio-angular models with slow rotational or translational diffusion promise unusual singular perturbation problems. Chapter 4 Spatio-angular order in populat ions of self-aligning objects 4.1 Introduct ion The topic of this chapter is the combined spatial and angular aspects of order that emerges in populations of interacting organisms. I generalize the three models introduced in Chapter 2 to the case of spatially inhomogeneous distributions. In order to do so, I first describe the random spatial motion of cells or molecules with additional linear terms. Second, to reflect the fact that the interaction between two elongated organisms depends on their distance apart as well as their relative angle, the interaction kernels are modified to include spatial dependence. The analysis in this chapter is restricted to linear stability analysis. This analysis permits separation of the leading modes breaking stability of the homogeneous state and gives a valuable hint about what kind of order may be expected. In the general spatio-angular case, the leading modes under certain conditions are spatially inhomogeneous. This is new in comparison with the previous two chapters. I concentrate on qualitative understanding of features of the interactions that are responsible for the appearance of the spatial or spatio-angular patterns. Theoretically, this is the case in the formation of spatially concentrated patches of an-gularly correlated cells (such as fibroblasts, for example), or of the growth of well-localized oriented filamentous network in cytoskeleton. The detailed mathematical description of 71 Chapter 4. Spatio-angular order in populations of self-aligning objects 72 this case is not possible. Nevertheless, the models investigated in this chapter give valu-able qualitative insights into the combined aggregation-alignment biological phenomena. 4.2 Descr ipt ion of the Mode l s Populations of objects are represented using a continuum description. Density distribu-tions that represent the population(s) of objects are functions of the space coordinate r £ D, orientational angle tt £ S of the object's axis, and time t. The spatial domain, D is either two dimensional (e.g., flat surface to which cells adhere in artificial in vitro growth conditions) or three dimensional (e.g., cells in vivo). In the two-dimensional case, the angle describing the direction of an object will be denoted 9 £ [—TT, 7r]. In the three-dimensional case, the angle in spherical coordinates is $7 = (0, 9) where 4> £ [0,7r] is the co-latitudinal angle, and 9 £ [0, 2ir\ is the longitudinal angle. The angular space is denoted by S in both two and three dimensions. Throughout this chapter, I consider only the cases in which the range of effective interaction between the objects is at least a few orders of magnitude smaller than the size of the domain. This is a natural consequence of the assumption about the contact-like character of interaction between the objects. But note that contact-like does not mean 8-like interaction, as they are commonly called in physics, since in that case, integral terms disappear from the models. Rather, the effective radius of interaction is of the same order as the length of the objects (usually much smaller than the size of the ensemble). This allows one to consider the spatial domain as infinite, and for the purposes of the linear stability analysis, to ignore boundary effects. Three separate models (the natural generalizations of the spatially-independent models introduced in Chapter 2) are studied. Chapter 4. Spatio-angular order in populations of self-aligning objects 73 4.2.1 Mode l I: Free and bound cells The dynamic behavior of the densities C and P is governed by the equations %(r,n,t) =[3CK*C + pPK*C-~(P, a*V ' ' (4.118) I §f(r ,f i ,<) = fi1AaC + fi2XC-(3CK*C-/3CK*P + -fP, where K(r — r', Q,£l') is the functional form of the rate with which an object at r',fl' turns to ft, and moves to r due to influence of object at r, ft; fi\ is the rate of random turning by free cells and ^2 is the rate of random walk (diffusion) of free cells. In these equations, AQ is the Laplacian in the angular variables. The operator A r is the Laplacian in spatial variables, i.e, ' £ + &, in 2D Ar 2 (4.119) A'* is a linear operator that depends on the orientation f), tt' and on the distance r — r' between the interacting objects. (K * f)(r,Sl,t) = [ dto' [ dr' K(r-r',n,n') f(r',n',t). (4.120) JS JD Details about the kernel, K(r — r',£l,£l') will be given below. The expression (4.120) is the influence of the distribution of material f(r', Q', t) on the angle 0 at position r. Thus a term of the form K * C represents the influence of the free object distribution at angle Q, and position r, and (3CK * C is the rate at which this influence results in free cells realigning and becoming bound cells at fl, r. The spatial Laplacian operator is used to capture a random motion of the objects which does not include pers istence in the direction of alignment, but rather pure ran-dom walk of the object's center of mass. In each experimental situation, there would be a specific type of stochastic process governing the motion of individual objects. Some Chapter 4. Spatio-angular order in populations of self-aligning objects 74 stochastic processes would lead to Brownian motion, while others result in a persistent motion in random direction. Further, random turning in general is described by more general operators, for example integral operators or their equivalents (see Murray (1989), section 9.5). Here we focus on the simplest situation, namely that of simple diffusion captured by the Laplacian operator. It will be seen below that there is a strong math-ematical reason for choosing the Laplacian, namely that its eigenfunctions are identical to those of the integral operator of the model. The form of the kernel K(r — r', ft, ft') in our model is of crucial importance and deserves special explanation. The dependence on the radius vector, r — r' between the centers of mass of the objects, follows from the fact that the spatial region is assumed to be homogeneous with respect to its influence on the objects. If this homogeneity is broken due to environmental bias or external force, then more general dependence K(r, r', ft, ft') arises. The convolution term describes an elementary process of alignment: two objects meet with initial directions ft, ft' and corresponding positions r,r'. With rate a, the objects continue moving with no interaction, and their directions are not changed. With equal rate (1 — <x)/2, the objects are attracted to either direction fl and position r , or to direction fl' and position r'. The factor (1 — a ) / 2 is absorbed into the constant /?. Consider the situation when the distance between the objects is L/2 (where L is the length of the object), and the objects are parallel to each other. Then we distinguish between two possibilities: (1) the directions of the objects are normal to the radius vector between their centers-of-mass, or (2) they are parallel to the radius vector. The value of K represents the probability amplitude of alignment. Experimental observations suggest that the greater the contact area between the objects, the greater is the probability of alignment. In the second situation, evidently this contact area is larger. For this Chapter 4. Spatlo-angular order in populations of self-aligning objects 75 reason, we can deduce that the value of K should be greater, though both values (r — r') and [Vt — Cl') are the same in the two cases. Thus, it is reasonable to assume that the interaction depends on the radius vector, and two angles of the objects' axes relative to the radius vector connecting them. For the sake of simplicity in this first investigation of the problem we will restrict attention to a somewhat less general form of the separable kernel, K(r - r', ft, ft') = K^Q - Q')K2(r - r'). (4.121) The strength of the interactions decreases as the distance between the objects grows. We assume that the effective range of interaction is of order L, where L is the length of the object. For the spatial part, it is natural to assume the form K2(r) = exp(^-), (4.122) where r is measured in units of L. For the angular part of the kernel we use the functions introduced in Chapter 2. Let us introduce the mass of objects (at all angles) at location r: M(r, t)= t (C(r, fl, t) + P(r, ft, t))dQ, (4.123) and the total mass of the whole ensemble, M= I M(r,t)dr. (4.124) J D It is easy to check that the total mass is conserved in the model. In the usual biolog-ical situation, the objects proliferate (e.g., cells undergo cell-division) so that the total mass increases. (This mass growth slows down and stops when the density of objects reaches some critical level.) As the objects proliferate slowly, we assume that M is the adiabatically varying variable. Chapter 4. Spatio-angular order in populations of self-aligning objects 76 For the purposes of analysis, it is convenient to recast the equations in a dimensionless form: f£( r , n, t) = CK * C + PK *C-aP, (4.125) f f (r, ft, t) = Cl AnC + e2ArC -CK*C-CK*P + aP. where a = j/(3 is the renormalized rate of exchange between bound and free cells, ei = f^i//3 is the renormalized rate of random turning by free cells, e2 — 1*21fl ls the renormalized rate of random walk by free cells. 4.2.2 Mode l II : Objects subject t o a l ignment and aggregation forces This model is more physically realistic, since it takes into consideration real forces causing real velocities. The motion of the object is taken to consist of the following two parts: (a) rotation about the center of mass with "angular velocity" ui = dfl/dt; (b) drift of the center of mass with linear velocity v = dr/dt. Velocities are assumed to be equal to corresponding forces: (4.126) v =Fr. The forces FQ, Fr are assumed to be conservative. The correponding potential functions are Fa =VnW, t ^ (4.127) Fr =VrV. The potentials W, V can be written in the form: W =W*C, (4.128) V =v*c, where W*, V* are the linear integral operators with the kernels Wi(£l)W2(r), Vi(fi)V2(r), respectively, where V2(r) = W2{r) = e x p ( - r 2 / 2 ) (4.129) Chapter 4. Spatio-angular order in populations of self-aligning objects 77 and Wi(il), Vi(Q) are the same as the kernel W^fi) from Chapter 2. These kernels represent the potentials created by a single object at (r', 0 ' ) acting on a single object at (r, fl). The interaction between the objects is effective at distances of order L between them and decreases quickly at large distances. It is possible to define fluxes (in both angle, J n , and physical space, J r ) as follows: Jfi = CUJ, (4.130) 3r =Cv. If we include both the divergence of the above convective fluxes and random motion, represented by the terms C^AQC and e2ArC for the spatial and angular diffusion, in a balance equation for the distribution, C ( 0 , r, t) we obtain the equation: BC — = e.AaC + e2ArC - V f i • J n - V r • J r . (4.131) Using the equations (4.126-131) in the above leads, finally to M o d e l II BC -^ = £ l A n C + e2ArC - Vn • (CVn(W * C)) - V r • (CVr(V * C)). (4.132) This equation describes the convectional drift of the objects in physical and angular space towards the points of the highest concentration, causing alignment and aggregation. 4.2 .3 Mode l III: Instantaneous aggregation and al ignment As in Chapter 2, the model again consists of a single equation. Mode l I I I : BC — = c.AuC + e2ArC. + C(Q(C) * C) (4.133) Here the integral term is defined as: Q(C) *C= I I L(C(r, SI) - C(r', fi'))G(r - r', ft - tt')C(r', Q', t)dSl'dr'. (4.134) Js JD Chapter 4. Spatio-angular order in populations of self-aligning objects 78 The kernel is now the product of two functions, one of which is responsible for the interaction dependence on the spatial and angular coordinates and has the same meaning and form as for the kernel K of the Model I of this chapter above: G(tt,r) = G1(n)e-r2/2. (4.135) The function, L(C(r,n) — C(r', 0 ' ) ) has the same meaning and properties as that in Chapter 2. 4.2.4 Eigenfunctions of the Operators As the spatial domain is infinite, the wavenumber q is a continuous variable. The complete set of orthonormal functions is {uq(r)zn(il)}, q£Rd, d = 2 , 3 , n = 0 , 1 , 2 , . . . , (4.136) where uq(r) ~ exp{i(q • r ) } , (4.137) and f ein9 in 2D, zn = { n = 0 , 1 , 2 , . . . (4.138) (Yn(<j>,e) in 3D. Yn are surface spherical harmonics (SSH). Our analysis is based on the following propo-sitions: 1. T h e set (4 .136-138) is the set of e igenfunctions of the operator A r wi th corresponding eigenvalues: 13, = -q2- (4-139) Chapter 4. Spatio-angular order in populations of self-aligning objects 79 2. The set (4 .136-138) is the set of eigenfunctions of the operator An wi th the corresponding eigenvalues: f - n 2 in 2D, an={ n = 0 , 1 , 2 , . . . (4.140) {-n(n + l) in 3D. Rigorous proof of the above two propositions can be found in any classical book on mathematical physics (concerning Proposition 2 see also Chapter 2). 3. T h e set (4 .136-138) is the set of e igenfunctions of the operators A'*, W*, V*, G*, with t h e corresponding eigenvalues: K(n,q), W(n,q), V(n,q), G(n,q), re-spectively. These eigenvalues are given by the following formula (in the case of K, for example): K(n,q) = kne-q2'2. (4.141) Similar formulae hold for W, V, G. The separability of the function A" is a direct con-sequence of the separable nature of the kernel K(il,r). (Here Kn are as in Chapter 2). Expression (4.141) follows from the fact that exp(—q2/2) is the Fourier transform of e x p ( - r 2 / 2 ) . 4.3 Linear stabi l i ty analysis In this section we consider only the 2D case. The 3D case is a clear generalization of our results. It is easy to see in all three models that the homogeneous solution (in both angle and space) is stationary. We are interested in the situation when pattern can arise. Our first step is to perform linear stability analysis to determine when the homogeneous solution becomes unstable. Mode l I Chapter 4. Spatio-angular order in populations of self-aligning objects 80 We consider the weakly perturbed homogeneous pattern: -P(r,n,t) where P, C obey the relations: — [PI c + [Pol Co uq(r) z„(0) e xt (4.142) M = P + C, P o + M ' M2 a + M ' (4.143) By substituting this into the system (4.125) and keeping terms linear in Po, Co, we obtain a condition under which A has a positive real part (implying growing perturbation, and thus, instability). We introduce the function A = Xn(q) which we call the Linear Growth R a t e , (LGR) . This function is obtained from the secular algebraic equation linking the values of A, n, and q, and represents the rate of growth of pattern close to the homogeneous steady state. Note that conservation of mass in all three models implies that the homogeneous solution is neutrally stable to perturbation by the mode n = 0, q = 0 (that is, Ao(0) = 0). The analysis is analogous to that in Chapter 2, and we get the instabi l i ty criterion: £2q' M o d e l II Substituting the pattern < - e m 2 + kne-fl\\ - kne-«2'2). (4.144) C(r,n,t) = C + C0uq(r)zn(n)eM, (4.145) into the equation (4.132) and keeping the terms linear in Co, we obtain the instabi l i ty criterion: e2q2 < - e m 2 + C(q2Vn + n2Wn)e^^2. (4.146) Chapter 4. Spatio-angular order in populations of self-aligning objects 81 Mode l III Substituting the pattern (4.145) into the equation (4.133), using the fact that at small C we have: L(C) = r]C + 0{C2), (4.147) and keeping the terms linear in Co, we get the instabil i ty criterion in the form: e2q2 < -em2 + nC2{\ - G ^ ' 2 ) . (4.148) The criteria for all three models lead qualitatively to the same conclusions. The most convenient way to investigate the conclusions of these inequalities is by graphical methods and by the comple te bifurcation diagram (see Figure 14), which represents the state of the system with all parameters fixed except for two governing parameters, ei and e2. These parameters are chosen because they are proportional to the angular and spatial diffusion rates, respectively, and we are interested in finding how the type of bifurcation depends on the relative values of these parameters. The shaded area in the parameter space is the region of stability of the homogeneous solution. The border of this region is the set of primary bifurcation points. The unshaded region is the region of instability. The dashed lines are the set of secondary bifurcation points. As the total mass of the system slowly grows, and t\, t2 are decreasing, we cross either the line labelled (A) or the line labelled (B). If we cross line A first then the mode corresponding to n ^ 0, q = 0 breaks the stability. Then the angular pattern evolves while the spatial distribution remains homogeneous, i.e., the objects orient but they do not aggregate. (See Figure 15.) The analysis of this situation is described in Chapter 2. We call this scenario A. If, on the other hand, we cross line B, then the modes n = 0,q « 0 (see Appendix C) are responsible for the instabilities so that the angularly disordered pattern with the spatial inhomogeneities Chapter 4. Spatio-angular order in populations of self-aligning objects 82 evolves. This means that the objects aggregate but they do not align. One would then see evolution of patches of objects within which the objects are not directionally ordered (See Figure 15.) We call this scenario B. The detailed calculations leading to these conclusions are given in the Appendix C. We now briefly describe a more general situation. Let us consider the kernel K(fl, r) such that its Fourier or Legendre coefficients have the form 2 2 Knexp(^-). (4.149) Here rn is the effective range of interaction (in units of V) within the angularly ordered phase. This means that spatio-angular dependence is more general and realistic than before. It can be shown that the kernel can be represented in the following form: K^^(2^- ( 4 ' 1 5 0 ) where the effective radius of interaction, K depends in a nontrivial way on the angle. Thus we control not only the strength of the interaction, but also its effective range as a function of angle. We note that not all such kernels with the Fourier or Legendre coefficients given above by expression (4.149) are positive. Since kernels in the model are transition probabilities restrictions must be imposed on acceptable values of Kn and rn. It can be shown (we omit details of the calculations) that the complete bifurcation diagram has qualitatively the same form as the one in Figure 14 but with displaced lines A and B. In particular, the bigger is rn, the larger the displacement from the origin of these lines. However, depending on Kn,rn, the diagram may change qualitatively in the vicinity of the intersection of lines A and B. (See Figure 16.) The new primary bifurcation (line C) appears corresponding to the mode n ^ 0, q ~ q^ / 0. This line in general corresponds to the largest value of rn. Chapter 4. Spatio-angular order in populations of self-aligning objects 83 If it happens that line C corresponds to the primary bifurcation or one of the secondary bifurcations, then the mode exp (in0) exp (iqr) in 2D and Yn(fi,) exp (iqr) in 3D is the leading one. Then spatio-angular order grows simultaneously. In this case, contrary to scenario B, the spatial density of the homogeneous state is not altered, but there exists long-range correlation between the axes of preferred orientation of the objects. In fact, the angle of preferred orientation changes periodically in space, with a characteristic period of order L, creating stripes (rolls) or spots (squares or hexagons). Schematically, the case of stripes is shown in Figure 15. This is strikingly similar to the experimental results of Tabony and Job (1990) in artificial polymerization of tubulin. We call this scenario C. 4.4 Discuss ion Mathematical results obtained in this section have the following qualitative interpreta-tion: we have two diffusion coefficients, the rotational one [X\ with dimension T~l and the translational one, u.2 with dimension l2/T. (I is the unit of length and T is the unit of time.) We can then introduce an "intensity of interaction" of order /3M (which is the same for both spatial and angular interaction and has the dimension T _ 1 ) and an "effective range of interaction" L. Three parameters with the dimension T _ 1 define time dynamics in this system: fj,i, H2IL2 •> (/3M). There will be an additional parameter, 7 in Model I. Then three parameters //17, ^27/ L2, and (/3M)2 determine the dynamics of the system. If IJ.2/L2 > (3M, over the effective range of the interaction translational diffu-sion destroys spatial order faster than the interaction between objects restores it. If simultaneously [i\ < /5M, the angular part of the interaction prevails over the rotational diffusion. In this situation scenario (A) is realized, and the objects align first, before any spatial order is evident. Otherwise, if Hi > /3M rotational diffusion destroys Chapter 4. Spatio-angular order in populations of self-aligning objects 84 angular order, and H2/L2 < /3M (spatial order is created), then we have scenario (B), i.e. spatial aggregation is evident before any kind of angular order is seen. For example, the usual situation for molecules is H1L2 ~ ^2, (see Landau and Lifshitz (1956)). In this case the detailed structure of the kernel determines the bifurcation scenario. However, if for one reason or another in a biological system the parameters, fii/L2 and ^2 are very different, then clear predictions can be made without knowing the detailed form of the kernel. The specific gaussian form of the kernels K, W, V, G is not crucial for the results of this paper. The only important thing is the symmetry of the angular dependence, and the fact that the spatial dependence is a decreasing function of distance which falls to zero on a length scale of the order of L, the size of the object. The smallness of the wavevectors responsible for the spatial bifurcation, implies a long length scale (relative to L) of the spatial pattern near criticality. This means that the length scale of the pattern may be comparable to the size of the spatial domain, and the effect of the boundaries could be significant. One of the implications is that irregular "droplets" containing large numbers of objects start to grow, and that a multitude of topological singularities can arise in the spatio-angular pattern (Gross and Hohenberg, 1993; Elsdale and Wasoff, 1976). This is similar to patterns formed in fibroblast cultures. It is of interest to consider more realistic forms of the kernels K(r — r', fi,fi '). Throughout the present chapter, I considered an isotropic domain. If there is anisotropy, the 0(n) symmetry is broken, and this must be reflected in the choice of kernel. I expect, further that nonlocal effects in the random motion of the objects will lead to nonlinear terms in the diffusion operators, which will cause important changes in the linear stability analysis (see, for example, Murray (1989)). In some systems in which there is substantial interest (for example, aggregations of animals such as fish schools, bird flocks, or herds of mammals) free motion of the Chapter 4. Spatio-angular order in populations of self-aligning objects 85 organisms is not limited to diffusion, but contains an ordered component of persistent motion, with some intrinsic velocity. (See Pfistner and Alt (1989)). Including this motion in the models would lead to hydrodynamic-like equations and would reveal a number of new phenomena. In all three models, we have at least two governing parameters proportional to the spatial and angular diffusion coefficients respectively. Due to the additional adiabatic growth of the total mass, the synchronous changes in these parameters causes one or an-other "growth protocol" (see Gross and Hohenberg (1993)) or a "developmental pathway" (Segel, 1984) namely, a sequence of the spatio-angular bifurcations. The first bifurcation is one of three possible types described above. The final pattern, however, may be formed by secondary bifurcations. (Moreover, the linear instability at q « 0 may not lead to any stability break. In this case, the mode responsible for one of the secondary bifurcations will be responsible for the patterns evolved.) One of the possibilities is that if spatial patches of objects are formed in an angularly disordered state, then the conditions for alignment are most favorable within the patches, and the axes of orientation of objects in the different patches are not correlated. If, on the other hand, patches start to grow in a partially aligned ensemble of objects, then patches themselves have elongated forms and the axes of orientation of the objects within the patches are correlated between neighboring patches. This phenomenon is seen in cellular automata simulations (EKE, 1990). It may happen that the spatial density of the homogeneous state is not altered, but that there exists long-range correlation between the axes of preferred orientation of the objects. In fact, the angle of preferred orientation would change periodically in space, creating stripes (rolls) or spots (squares or hexagons). This kind of order was observed experimentally (Tabony and Job, 1990) in condensation of tubul in in vitro (tubulin is Chapter 4. Spatio-angular order in populations of self-aligning objects 86 one of three most abundant polymers in the cytoskeleton). It is interesting to note that there are two groups whose work on liquid crystals is vaguely in the flavour of our approach. The first due to Chandrasekhar et al. (1970), is a Monte-Carlo simulation of the finite system of elongated molecules represented by objects (ellipses) which undergo thermal transitions restricted by excluded volumes (the ellipses are not allowed to overlap). It is found that a sequence of two transitions - po-sitional and orientational - is observed as the density of the objects increases. For less elongated ellipses, when the homogeneous steady state is disturbed, first the positional transition occurs (patches of spatially ordered ellipses formed) and only then the rota-tional transition. For more elongated ellipses, the order of the transitions is reversed. In these simulations, direct interactions between molecules is absent. However, due to the effect of excluded volumes it is most geometrically favorable for them to take on some degree of spatio-angular order. There is clear analogy of these results on the sequence of phase transitions with our scenarios (A) and (C). An example of an approach that resembles our model in the field of liquid crystals is the paper by Greco and Marrucci (1992). They also postulate a dynamic system of nonlinear PDE's to describe the inter-actions of rod-like molecules, and their distribution over space and angle. The details of the equations, and the techniques and results of the analysis is, however, quite different from our own. Cellular automata modelling of a similar system was described in EKE (1990); Er-mentrout and Edelstein-Keshet (1990). In these models, persistent motion of the objects at random directions, rather than spatial Brownian motion was considered. Furthermore, hard-core repulsion between the objects was introduced. Despite these differences, the qualitative picture of the patterns evolved bear similarities with the qualitative predic-tions we make in this paper. The range of interaction is virtually greatest in the aligned Chapter 4. Spatio-angular order in populations of self-aligning objects 87 state, so that scenario (B) is not realized. In some cases, it was clearly seen that patches start to form in the aligned phase, resembling scenario (A). The form of the patches is elongated in the same direction as the axis of orientation of the objects within the patches. In other cases (when the density of the objects was higher), patches are formed from aligned objects, but different patches have different orientations. This suggests that scenario (C) can be realized. We emphasize that cellular automata modelling produces pictures of well-developed patterns, while our paper here deals with the initial phase of this process. No direct correspondence can be made for this reason. Chapter 5 Partial integro-differential mode l for aggregation in a bacterial swarm 5.1 Introduct ion 5.1.1 Swarming behavior of Myxobacter ia Myxobac ter ia are procaryotes with a multicellular life cycle (Shimketz, 1990). They are distributed worldwide due to the resistance of the myxospores to desiccation and freezing. They are a common component of the soil. Virtually every pinch of soil all over the world contains at least a few bacteria of this type. Myxobacteria provide a vast array of antibiotics with unusual chemical structures, some of which may prove to be clinically useful. In this chapter, the collective behavior of Myxobacteria, their swarming, their rate of group motion, and patterns (such as ripples) formed in their colonies will be considered. A description of the properties of an individual, and their mutual interaction is essential. I want to demonstrate how patterns observed in a Myxobacterial swarm could be understood from the level of interactions between individuals. An integro-PDE model will be used to investigate rippling. The mathematical goal is to explore analytically the rich variety of patterns described by an integro-PDE. The biological goal is to explain which specific details of chemotaxis and direct interaction are responsible for pattern formation. A bacterium has the shape of a semi-rigid rod, bending just slightly, with average length 5 fj,m and 0.5 //m wide. The cell lacks motility organelles (such as flagella, for instance), and performs so-called gliding motility. When in motion, the bacteria glide 88 Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 89 back and forth, changing their direction often. Though the instantaneous velocity of gliding is significant, the average velocity is of the order of 2-3 /xm/min due to frequent turning. The angular orientation of a single cell undergoes a slow change with a rate of about 3° per minute (Kaiser and Crosby, 1983). The mechanism of gliding motility remains elusive. Myxobacteria move in a dense extracellular matrix, composed of pro-tein and carbohydrate (Behmlander and Dworkin, 1994). The cells spend a large part of their metabolic energy on creating and maintaining this extracellular matrix (often inappropriately referred to as slime). Gliding cells often move in large, organized groups, known as swarms, in which bac-teria maintain continuous contact with their nearest neighbors (Dworkin and Kaiser, 1985; Reihenbach and Dworkin, 1981; Reihenbach, 1984, 1986; Shimketz and Kaiser, 1982; Kaiser and Crosby, 1983; Shimketz, 1990). The spreading of the swarm results from cell motility, not cell division or growth. Extensive research has revealed that two separate genetic motility systems control the cell movement (Shimketz, 1990). System A controls the so-called "adventurous" motility, enabling an individual bacterium to glide on its own. In close proximity to another cell, the system S (for "social" motility) is activated, governing cooperative movement. Social behavior of Myxobacteria depends on three interrelated properties: - gliding motility, - secretion a variety of substances, - cell-cell contact mediated interaction. The reasons for swarming behavior in bacteria are still unclear. One hypothesis is that it is advantageous to the cells for the purposes of motility or feeding to share the production of secreted chemicals (such as enzymes). This can be called a "wolfpack" phenomenon in which a group is more effective than an isolated individual. Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 90 The average velocity of a swarm depends on the bacterial density (Kaiser and Crosby, 1983). (A typical value of the density is 106 cells/cm2.) The maximal spreading rate is achieved for a dense swarm in which the center-to-center distance between cells is of the order of 5 /mi. (This rate stays approximately the same for denser swarms.) At low densities the spreading rate is several times smaller and grows linearly with density. The formula, suggested in Kaiser and Crosby (1983) for v, the local velocity of a swarm, v = V(l- e-c'c), (5.151) approximates the experimental data well. Here V is some parameter - limiting velocity, C is the density of a swarm, and c is some parameter representing the density at which the linear growth of velocity slows down. The bacteria are aligned within a swarm; their average rate of reorientation is even smaller than that of an individual cell. The behavior on the edge of a swarm is very peculiar - cells initially advance rapidly at the edge, and then, as if pulled back by an invisible thread, reverse and rejoin the swarm. Sometimes the edge consists of "flame-like" sharp protrusions of the most active, "adventurous" bacteria. In other cases the average density at the front of the propagating swarm decreases smoothly. There is experimental evidence for the existence of both short-range and long-range interaction in Myxobacteria. When two single cells meet, they perform a distinct ritual; they "rub" against each other, and then abruptly change their initial directions of motion (Reihenbach and Dworkin, 1981). Bacteria are also thought to communicate by chemical transporters via fibers of the extracellular matrix (Behmlander and Dworkin, 1994). These fibers extend for a few cell lengths, thus increasing the spatial range of interaction. A one-dimensional model of a Myxobacteria swarm considered by Pfistner (1989) deals with densities of bacteria of two sorts - those gliding to the right, and those gliding to the left. The speed of gliding was taken to be constant and the interchange between Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 91 the two groups due to random turning and non-local interaction was assured. This model consisted of a system of two integro-PDEs, and has yet to be treated analytically. It was shown numerically that these complex processes provide the existence of a travelling band. On a more realistic level it would be desirable to investigate a two-dimensional model in which individual bacteria can adjust the direction as well as the speed of movement. A sketch of a model of this type was developed by Pfistner and Alt (1989). Even the linear stability analysis of this model would be non-trivial. 5.1.2 Rippl ing The phenomenon of rippling was discovered by Reihenbach in the 1960's in the course of time-lapse microphotographic studies. Myxobacteria are the only procaryotes per-forming this peculiar rhythmic behavior. The rippling swarm shows radial veins of non-homogeneous density with trains of travelling waves - periodically spaced ripples (Dworkin and Kaiser, 1985; Shimketz, 1990; Shimketz and Kaiser, 1982). Cells are con-centrated at the ridges and depleted from the troughs. In rippling, the average density of cells is as high as 4 • 107 cells/cm2. The swarm is about 5 fj,m thick, and thus the volume density of bacteria is then 8 • 1010 cells/cm3. Cells in the crests are more highly aligned than those in the troughs. An intriguing observation is that their axes are not normal to the ridge, but subtend an angle of approximately 40° to the ridge direction. The distance between ridges and the width of the ridges is of the order of 50-100 fj,m. The ridges move synchronously in the same direction in parallel tracks, 100 to 200 jiva wide. The velocity of the ripples is 1-3 / im/min, close to that of a single bacterium. Individual ripples exist for 15-20 min, though the field of ripples as a whole persists for up to 5 hours. Due to their skewed orientation individual bacteria gradually drift to the lateral side of the track and cross from one track to another. There is little net displacement of the bacteria in Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 92 spite of the persistent propagation of the ripples. Both A and S motility systems are essential for rippling. Mutants that are deficient in one or both, do not exhibit this rhythmic behavior. The causes and role of rippling in the Myxobacteria life cycle is not yet understood, but is a phenomenon known to be prevalent in fruiting body formation. Rippling is often observed in situations of substantial preying. One of the genes required for rippling is also essential for development (Shimketz, 1990). Also, an opinion exists that this type of aggregation is perhaps some artificial degenerate swarming or a phenomenon similar to that of Dictyostelium aggregation (Shimketz, 1990). 5.2 Non- loca l aggregation mode l . 5.2.1 I D mode l The basic model I suggest is closely related to Model II of Chapter 4, a model, in which objects, in this case bacteria, influence one another's motion through direct or indirect effects. In a ID situation it is possible to derive the terms of the final model from a deeper consideration of microscopic effects, for example from some (hypothetical) underlying chemistry. I start with the continuous Eulerian approach. In ID we introduce the linear density of the swarm C(x, t) (cells/unit of length) as a function of the coordinate x and time t. I assume that all cells align fully and glide along the x axis. One way of deriving a model is to assume that bacteria secrete two types of diffusible chemicals, chemoattractant and chemorepellent. A second way in which some "direct" interactions are assumed, instead, will be explored in the next section. Let the densities of the chemoattractant and chemorepellent be si(x,t) and s2(x,t), respectively. The spatio-temporal dynamics of these three densities would then be given by the following PDE Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 93 system: dt dx2 1dx\dx J 2dx\dx / ' . . dsi d2Si t i ~ = di-^-j- + ciiC - bsi, i = 1,2, e; < 1. The first term in the RHS of the first equation describes cell diffusion with diffusion coefficient Di. The next two terms represent drift up the gradient of the chemoattractant and down the gradient of the chemorepellent, respectively. {Ax and A^ are parameters for the magnitude of the chemotactic response.) In the second equation the e,'s are relaxation time constants of the chemicals. I assume that the chemicals have very fast dynamics, and that their distribution is quasi-stationary, following the bacterial density distribution. This means that et's are small parameters. Chemical diffusion with diffusion coefficients d\ is described by the first term. The second term denotes the secretion of the chemicals by bacteria with corresponding rates a;'s. The last term represents decay of the chemical substances with rates b{. This model is a simple generalization of the earliest slime mold aggregation model by Keller and Segel (1971). In the so-called adiabatic approximation, we set the LHS of the second and third equations of the system (5.152) equal to zero. The quasi-stationary spatial distributions of chemoattractant and chemorepellent are given by: & - %* = -%C. (5.153) This linear equation is easily solved by the Green's functions method (Zwillinger, 1992). The solution reads: *i = V%uij_00 e x p ( _ v f lx - *'l) C(x')dx'. (5.154) The Green's functions method allows one to express the chemical distribution in the form of a convolution integral. The integral kernel here represents the spatial spread of Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 94 a chemical emitted by a cell at the point x' to the point x. Note that this is exponential decay over distance. Consider the new set of parameters: *=*&> B = ^ a = ^ b = ^ (5-155) and the function K(x) = sign(x) [Ae~alxl - Be~b^) . (5.156) (The function sign is equal to + 1 when x is positive and to —1 if £ is negative). Note that K is an odd function. The parameters (5.155) have the following meanings: A = amplitude of the repulsive part of the interaction, B = amplitude of the attractive part of the interaction, a~x = the spatial range of the repulsion, b~l = the spatial range of the attraction. Plugging the solutions (5.154) into the first equation of the system (5.152) and using (5.155-156), we obtain the following equation for cell density: ID Model dC d2C d r+°° -5T = A T T T ~ ir(vC)i V = (K *C)= K(x - x')C{x')dx\ Ot UX OX J — OO /~ * rry\ (5.157) K(x) = sign(x) (Ae~aW - Be~hW) . This model is similar to Model II of Chapter 4. This single diffusion-advection equa-tion governs swarm dynamics in the limit of fast chemical dynamics. As the drift velocity, v, is density dependent, the equation is non-linear. Moreover, the expression for the ve-locity is given by a convolution integral term and is thus non-local: the direction and magnitude of the drift at a point x depends on swarm densities at all other points x'. This dependence is expressed by the convolution kernel K(x) (further referred to as the kernel ). Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 95 The expression for the velocity suggests the following alternate interpretation: Imag-ine that a cell at the point x' affects the velocity of a cell at x by "an amount" K(x — x'), without any intermediate substances. Assuming linear superposition of the influences of all cells we get the expression for the velocity of a given cell. This implies that direct interaction can be described within the same mathematical framework. Further, this suggests the following helpful mechanical analogy. It can be assumed that bacterial ve-locities are proportional to some fictitious forces. This assumption is common in biology at low Reynolds numbers, see Mogilner and Edelstein-Keshet (1995a). Then K(x — x') is proportional to the force with which a cell at the point x' acts on a given cell at the point x. Without loss of generality, we can take the proportionality coefficient as unity; then K(x — x') is the corresponding force. Note that the function K{x) is odd, so that in each pair-wise interaction, the distance between a pair of cells changes, while the center-of-mass of the pair remains intact. For the whole swarm this leads to conservation of the center-of-mass of the swarm (see Landau and Lifshitz (1976) for the mechanical analogy). Biologically this means that without external bias (not yet included in the model) there can be no macroscopic swarm motion. The first term in the expression (5.156) can be treated as a repulsion and the second term as an attraction between a pair of cells separated by a distance x. A and B are the corresponding strengths of the interactions, while a~l and b~l are the spatial ranges of the interactions. 5.2.2 2 D Mode l In 2D, a full derivation of the model from some underlying (anisotropic) chemotactic system is problematic. However, the direct-interaction system of Model II of Chapter 4 can be applied as is, by simple generalization to 2D. We introduce the swarm density Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 96 C(x,y,t) as a function of x and y coordinates and time t. All cells are assumed to be densely packed and thus completely aligned in the x direction. I also assume gliding motion is predominantly along the x axis. Then the following equation governs the 2D density distribution of the bacterial swarm: Here the first term describes drift along the x axis with average velocity v. Random motion in the x direction, the stochastic part of the gliding system, with a corresponding coefficient D\ is given by the second term. The third term describes slow random motion in the y direction with coefficient D2- This "diffusion" is observed experimentally and results from stochastic exchange of positions of neighboring cells. The rod-like shape of the bacteria and the dense cell packing imply that diffusion in the y direction is much slower than that in the x direction (D\ » D2). The expression for the drift velocity is given by a convolution integral term: dx' / dy'R(x -x',y- y')C(x', y'). (5.159) -00 J—00 The kernel is now a function of the coordinates (x,y) and (x',y') and represents the contribution to the velocity of a cell at (x, y) by that at (x', y') (Figure 20). This influence does not depend simply on intercellular distance because the swarm is anisotropic. We chose the simplest factorized form for the kernels: R(x -x',y- y') = (K(x - x') + sKe(x - x'))W(y - y% (5.160) where s is some dimensionless parameter, and the form and meaning of the kernels K(x) and Ke{x) are described below. The dependence on the distances x — x' and y — y', rather than on the more general forms (x,x') and (y,y'), are chosen to reflect translational invariance of the system in each of the directions x and y separately. Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 97 (K(x) + sKe(x)) reflects dependence of the interaction on the separation along the x axis (axis of motion of the swarm), and W(y) similarly represents dependence on the y separation. The kernel (K + sKe) governs the motility in a biased case. It consists of two terms, the first an odd function for the non-biased case and the second Ke(x) (an even function) is not responsible for convergence or divergence of the cells, but for an average net drift. Ke(x) describes the motion of the swarm as a whole for example, when some external bias like an uneven distribution of nutrients causes directed motion in the x direction. This follows from the fact that the convolution with Ke(x) is non-zero for constant swarm density. Note that here we do not consider the density-independent com-ponent of the velocity. The corresponding term can always be omitted by transformating to a coordinate system moving with the given velocity. The parameter s represents the relative strength of biased and non-biased interactions. The non-biased case corresponds to the value s — 0. For the function K(x), we use the ID expression (5.156). We do not specify the expression for the function Ke{x) other than stating that it is some even, non-negative function with finite support. The kernel W(y — y') is the factor weighting the strength of interaction between a pair of cells as a function of the distance in the y direction. By symmetry of left and right directions along the y axis, the function W(y) should be even. Furthermore, it is reasonable to assume that the strength of interaction decreases monotonically with increasing distance, \y — y'\. Two functions having these properties will be considered: W1(y) = Jre-y2'\ W2(y) b M < i (5.161) 0, |y| > 1. In the first case, interaction gradually decreases in the y direction, while in the second, interaction is constant for small distances, and then abruptly disappears for distances Chapter 5. Partial integro-difFerential model for aggregation in a bacterial swarm 98 greater than one "unit". I will call these cases "continuous" and "finite-range" interac-tions respectively. In summary, the two-dimensional model has the form shown below: 2 D Mode l dc d . „. n d2c „ d2c dt dx dx2 dy2 ' »(*> v) = / - £ dx' / -« W(K{x - x>) + sKe(x - x'))wh2{y - y')C(x>, y>), K(x) = sign(x) (Ae~a\x\ - Be'^) , ( 5 J 6 2 ) h \V\<1 W^y) = -fee-*l\ W2(y) = 0, | y | > l As shown in the next three sections, the Fourier transforms of the kernels play a decisive role in pattern formation. Consider, specifically, the transforms: f+oo A r+oo Ke(q) = I Ke(x)cosqx dx, K(q) = / K(x) s'mqx dx, (5.163) J —CO J —CO r+oo W(q)= W(x) cos qxdx. (5.164) J —CO For the functions K(x), W\(y), W2(y) (5.162) used in the present chapter, we have: W,(q) = e-*'l\ WM = —- (5.166) The complete 2D picture is as follows: There are two regions surrounding each bac-terium: The given cell is attracted (repelled) by all cells in region 1 and repelled (at-tracted) by cells in region 2 (see Figure 20). (Which situation is realized depends on the Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 99 model parameters.) Interactions with cells outside these regions are exponentially small. Also, the average net drift velocity of a cell depends on a weighted local swarm density. This single integro-PDE model (equation (5.162)) leads to a great variety of patterns, as descibed in the next 3 sections. 5.3 Travelling wave solution. I investigate predictions of the model under the assumption of slow spatial variation of the density (meaning that the density distribution does not change significantly over the spatial range of the kernel). The standard procedure in this case is to expand the function C(x,y,t) in a Taylor series and to keep the first few terms (Murray, 1989; Lewis, 1994). This approach is presumably valid for a rough description of the large-scale behavior. The phenomenon of the Myxobacteria swarm suggests looking for the solution in the form of a single travelling plane wave propagating in the x direction. Thus, we look for a solution of the 2D equation in a form independent of y: C — C(x,t). Expanding the swarm density function in a Taylor series: C(x', t) = C(x, t) + (x' - x ) ^ ^ - + ..., and keeping the first two terms, we obtain the following expression for the local velocity from equation (5.162): dC v^s^C-Vi — , (5.167) ox K(x)x dx = l imK(q)±, (5.168) -oo q-+0 where the expressions Ke(0) and K(q) are given by the definitions (5.165). The coefficient Vi represents the density-dependent part of the biased local velocity of the swarm. The coefficient V2 characterizes the amplitude of the local swarming effect. Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 100 Note that this approach turns a non-local expression into a local one. Plugging (5.167) into the equation (5.162), we get the reduced PDE: (The z/-derivative disappeared because of the plane-wave character of the solution.) A travelling wave solution would have the form: C{x,y,t) = U(x-Vt) = U{z), where V is the travelling wave velocity which is to be found. Plugging this form into the equation (5.169) and integrating once leads to the ODE: -VU + (sV1U-V2U')U-D1U'= H, (5.170) where H is an integration constant and ' denotes a derivative with respect to z. Let C_ and C+ denote densities far ahead and behind the swarm front, i.e., such that l im^-oo U = C- and lim2_^oo U = C+ = 0 (assuming that there are no bacteria far ahead of the swarm). Using these limits, one finds that H = 0 provided the velocity of propagation of the swarm is V = ViC-. Then ^= 5 V i^ W ' (5-m) ^-InU + (Dl + V2°^ ln(g_ -U) = h- sKz, (5.172) where h is a constant of integration. (The value of h is unimportant: changing h just shifts the travelling wave solution as a whole along the x axis.) The asymptotes of the travelling wave have the form: z» nftr> u - h+ e x p ( SD\ " Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 101 where h+,h- are constants. This travelling wave solution describes the propagation of a swarm of persistent shape, resembling that observed (see Figure 21). The velocity of the front edge of the swarm V is proportional to the bacterial density in the center of the swarm (V = ViC_), which agrees with observation at intermediate densities (see the phenomenological expression (5.151)). If the individual velocity does not depend on the density (s = 0), the non-trivial travelling wave solution is absent in this approximation. The solutions (5.173) are good approximations only in the case — — < mtn(a, b), < min(a, b) L>\ L)\ -f- V2^y-(i.e., the front of the travelling wave changes slowly on the length scale of the range of interactions. Recall that the parameters a _ 1 , 6 _ 1 defined in (5.155) represent the ranges of repulsive and attractive interactions respectively). It follows from these inequalities that when diffusion becomes too slow, or when the cell density grows too high, the local approximation is no longer valid. The reduction of the integro-PDE to the PDE is also invalid if V2 < 0 and D\ + V2C- —> 0. Biologically, when the local swarming effect becomes strong enough, dispersal is completely suppressed and the swarm becomes too dense. Mathematically, the problem then becomes ill-posed (see Alt (1985)). It will be shown that the last condition leads to a bifurcation to a dense bacterial patch which cannot be described within the framework of the local approximation. Finally, the existence of periodic solutions of the original model will be demonstrated later. The period of these patterns is of the order min(a~1,b~1). In this case, all terms in the Taylor series approximation are of the same order so that even if we keep higher terms of the series, results are bound to be questionable (see also Grunbaum and Okubo (1994)). In summary, the reduction of the integro-PDE to the PDE successfully gives a travelling wave solution when diffusion is sufficiently fast and swarm density is low. Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 102 5.4 I D Linear and non-linear stabil i ty analysis. 5.4.1 Linear stabi l i ty analysis I now consider the problem (ID Model (5.157)) on an infinite domain, and investigate behavior far from the front edge of the swarm and far from any natural boundaries of the experimental setting. The initial swarm density is assumed to be constant throughout an infinite domain. Linear stability analysis of the ID model is used to examine possible patterns. Let C be the homogeneous stationary solution of the equation (5.157). Consider small perturbations of the homogeneous solution C = C + c. Plugging this expression for C into the ID model equation (5.157) and keeping only terms linear in c one can easily derive the linearized equation: dc d2c - d (Note that (K * C) — 0 by symmetry properties of the kernel K(x)). Harmonic plane waves (of the form c = eXtetqx, where q is a continuous wave vector, and A is the cor-responding linear growth rate) are eigenfunctions of both the differential and the linear integral operators on the RHS of equation (5.174). Plugging this form into (5.174) we obtain the dispersion relation: A = -Cqk(q) - Dxq\ (5.175) where K{q) is given by the formula (5.165). The second term on the RHS of the dis-persion relation (5.175) is negative. Its absolute value increases as the value of the wave vector grows. Therefore, the homogeneous steady state is stable to short-range excita-tions. On the other hand, the homogeneous state is neutrally stable to very long-range perturbations by the translational invariance of the model. Thus, one may expect insta-bility to take place when the first "interaction" term in the RHS of (5.175) is positive Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 103 and large enough at a finite wave vector q = qc. It is convenient to introduce the bifurcation function \{q): % ) = o — = 2 a 2 - T ^ — o ' (5.176) 2# gJ + b2 q2 + cr Note that this function is even. Its properties are described in Appendix D and summa-rized in Table II. The linear growth rate has the form: A = ( 2 ( 7 % ) - A ) q2. (5.177) The instability occurs in the case: ReX>0, % ) > ^ - (5-178) (Note, that the bifurcation function does not depend on the diffusion coefficient). The most convenient way to investigate (5.170) is by a graphical method. To investigate the character of possible bifurcations, we fix all parameters except D\,B,b, which are the governing parameters. Thus, we fix the strength of the repulsion A and its range a - 1 , and vary the strength of the attraction B and its range 6 - 1 . It is convenient to introduce as the governing parameter space the quadrant B > 0, b > 0, and to subdivide it into 10 sectors with the 5 curves shown in Figure 22. In different sectors the interaction kernel (5.156) and the bifurcation function (5.176) have qualitatively different forms, which lead to different possible patterns. The curves are defined by the equations: (i) : b = a; (ii) : b = a(BfA)1!4; (Hi) : b = a(BfA)1'2; (5.179) (iv) : b = a(B/A); (v) : B = A. (The reason for so many sectors is that the qualitative forms of the interaction kernel and the bifurcation function are very sensitive to the values of the governing parameters. It is possible, in principle, to reduce the number of the sectors, but then the analytical Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 104 description would become cumbersome). In what follows, we analyse the conclusions from the instability criterion (5.178) at values of the governing parameters in each of the regions. If the values of the governing parameters B and b are such that the bifurcation function is not positive, then the instability criterion (5.178) cannot be valid, and the homogeneous steady state is stable. If, on the other hand, the bifurcation function is positive at some values of the wave vector q, then an instability can occur. The leading unstable mode is characterized by a value q = qc at which the bifurcation function has a global positive maximum. (Since the bifurcation function is even, bifurcation always occurs symmetrically at ±g c . Below we will drop the sign ± with the understanding that this symmetry exists.) The necessary condition for a stability break at this value of the wave vector is that the diffusion coefficient becomes small enough, that the RHS of the instability criterion (5.178) is smaller than the value of the bifurcation function at the critical value of the wave vector. The bifurcation function is plotted at positive values of the argument for different val-ues of the governing parameters B and b in Figure 23. In order to connect the qualitative behavior of the system with the properties of inter-bacterial interaction, the interaction kernel is also plotted in Figure 23 (positive values of the interaction kernel correspond to repulsion, and negative ones to attraction). Now I can make the following statements about the stability of the homogeneous state based on the qualitative analysis of the bifurcation function: Region 1, 2 and 10 (B < A,b > a(B/A)1/2): For model parameters from regions 1 and 2 there is effective repulsion at short distances and attraction at larger distances. (Quantitatively, repulsion is much larger than attraction.) Parameters from region 10 Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 105 correspond to pure repulsion. The bifurcation function is negative everywhere, the in-equality (5.178) cannot be valid, and the homogeneous steady state is stable. Regions 3, 4, 5 and 6 (B < A,b < a(B/A)1/2 and B > A,b < a{BjA)1lA): Parameters from regions 3 and 4 correspond to effective repulsion at short distances and attraction at larger distances. Region 5 represents pure attraction, and that from region 6 to attraction at short distances and repulsion at larger distances. Quantitatively, in all these cases, attraction is much stronger than repulsion. The bifurcation function has a positive maximum at the origin. As the diffusion coefficient decreases, break of stability always occurs in a small vicinity of the critical wave vector qc. Thus, in this region, at some critical diffusion coefficient, long-wave harmonics break stability of the homogeneous state (qc ~ 0). Growth of pattern on a large spatial scale can be expected. Observe that in this case the instability criterion can be written in the form D\ < —V2C, for V2 as in (5.168). The sign of V2 is the sign of the response of a cell to the local gradient of swarm density. (For example, V2 negative and large means that the cell is drifts up the gradient.) From the present analysis, it follows that in this case we can expect some non-periodic large-scale pattern. However, as shown in section 3, the local approximation then fails. Regions 7, 8 and 9 (B > A,b> aiyBjA)1'A): For parameter values from these regions, interaction is effective attraction at short distances and repulsion at larger distances. Quantitatively, repulsion is stronger than attraction. The bifurcation function has a positive maximum at qc = (B^a2 - Ah2)' (A* - B*)~* . (5.180) (The extremum of the bifurcation function away from the origin, if it exists, is always given by this expression, see Appendix D.) As the diffusion slows down, the break of stability occurs at this qc ^ 0, and we expect the appearance of periodic patterns. Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 106 The results of linear stability analysis are summarized in Table 2. The biological conclusions from this analysis are as follows: 1. If effective repulsion is very strong, the homogeneous steady state is stable. 2. If effective attraction is sufficiently strong, bifurcation from the homogeneous state caused by long-wave harmonics occurs in the limit of slow diffusion or high swarm density. 3. If there is a balance between attraction and repulsion (such that attraction is short-range and repulsion is long-range), periodic pattern may emerge. To establish whether the patterns actually grow, we have to perform weakly non-linear stability analysis. This is discussed in the next section in the case of periodic pattern. Unfortunately, in the situation that the critical wave number qc ~ 0, the non-linear analysis becomes very non-trivial (Gross and Hoenberg, 1993; I. Aronson, private communication), and will be performed elsewhere. It is worth mentioning that one of the conditions for periodic patterns in a Turing-type system is the presence of long-range inhibition and short-range activation (Segel, 1984). This may be analogous to our model. 5.4.2 I D Non-l inear analysis In this section I apply the synergetics approach to weakly non-linear stability analysis of ID Model (5.157). The treatment is similar to the one already described in Chapter 2. Amplitudes of the growing modes are assumed to be independent of the spatial coordinates. Instead of investigating all unstable modes characterized by wave vectors from a narrow band (the so-called "small band excitations") only one "fastest" mode is considered. In this way, slow spatial modulations of the periodic pattern are neglected. (See the detailed discussion of this approximation in Haken (1977); Gross and Hohenberg (1993).) Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 107 We expand the density distribution in the form: C = l + Gcos(qcx)+ ^2 g(q) cos (qx), (5.181) where G is the amplitude of the unstable mode, qc ^ 0 is the critical wave number and g(q) are amplitudes of the stable modes. Without loss of generality, take (7 = 1, and use the sum rather than the integral for summation over the wave vectors. (Formally this corresponds to a periodic spatial domain in the limiting case of a large period approaching infinity.) Then the space of wave vectors is discrete, and the distance between adjacent vectors diminishes to zero. (See the discussion in Haken (1977).) The stable amplitudes obey the equation: 9(q) = \(q)9(q) - I £ K(qW) 9(qW) 9(q{2)) S(q ~ q" - q^). (5.182) Z g ( l ) , g ( 2 ) (This equation is obtained by plugging the Fourier integral for the density distribution into equation (5.157) and gathering all terms corresponding to the same wave vector. Here g(1) and q{2) are pairs of wave vectors corresponding to interacting modes, and X(q) is the linear growth (or rather relaxation) rate (5.175).) By the "slaving principle" we set the LHS of equation (5.182) equal to zero, and then assume that G ^> g(q)- This allows us to expand the stable amplitudes in an asymptotic series in powers of G. The stable stationary amplitude g(2qc) of the wave corresponding to 2qc, is of order G2. It is given by the asymptotic expression All other stable amplitudes are of order G3 and can be neglected in the first approxi-mation. Using this last expression and (5.182), we obtain the generalized Ginsburg-Landau ODE for the unstable amplitude: 2 G = SXG -pG3, P= -r£-rK(qc) (k(qc) - K(2qcj) > 0, (5.183) Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 108 where SX — X(qc) is a small positive bifurcation parameter. These expressions give stationary amplitudes of orders: G-V8X, g{2qc)~8\, g{q)~8\3/\ q ^ qc,2qc. (5.184) It follows from expression (5.181) that close to the bifurcation point, there is a small stable stationary amplitude of the nonhomogeneous part of the solution. Thus, the weakly non-linear ID stability analysis demonstrates that when linear analysis predicts periodic pattern, this pattern, indeed, emerges close to bifurcation. Let us note, however, that the sign of the Landau constant p is of crucial importance. The analysis is valid for positive values of p only, a condition that is satisfied for our choice of kernel. It may happen, though, that for some other reasonable choice of kernel, p < 0. In that case the non-linear analysis presented here would be invalid. The problem of the bifurcation scenario in the latter case remains open. The possibilities in this situation include: 1. Higher order non-linearities limit growth of the unstable mode, and, still, the periodic pattern appears. 2. Non-stationary, oscillating, or chaotic patterns grow (D. Wollkind, private com-munication). 5.5 2 D Linear and non-l inear stabil i ty analysis. 5.5.1 Linear stabil i ty analysis I now consider the 2D Model (5.162) on an infinite plane. As in the ID case, the homogeneous state is stationary. We are looking for the perturbed homogeneous state in the form C = C + c. Plugging this expression into model equation (5.162) and keeping only terms linear in c one can easily derive the linearized equation: Be - 3 ~ dc d2c d2c | = _C | (H.c)-(fl.C)| + B l ^ + f l , W (5.185) Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 109 Harmonic plane waves (of the form c = eAte ,(?i :c+,2S '), where q = (qi,q2) is a continuous wave vector with x component q\ and y component q2, and A is the corresponding linear growth rate) are eigenfunctions of both the differential and integral linear operators on the RHS of equation (5.185). Plugging this form into (5.185), we obtain the dispersion relation: A = -isCqi(ke{0) + Ke(qi)W(q2)) - Cqik(qi)W(q2) - Diq\ - D2ql (5.186) where the Fourier transforms of the interaction kernels are given in (5.165-166). The last two terms on the RHS of the dispersion relation (5.186) are negative. Their absolute values increase as the modulus of the wave vector grows. Therefore, the 2D homogeneous steady state is stable toward short-range excitations. As in the ID case (see Section 5.4.1), the homogeneous steady state is neutrally stable to very long-range perturbations. Instability is expected, then, at some finite value of q = (qic, q2C) when the second (real) "interaction" term on the RHS of (5.186) is positive and sufficiently large. The meaning of the imaginary term in the dispersion relation becomes clear if we write the expression for the profile of the plane wave excitation in the form: c _ e^tei{qix+q2y) = exp[(~Cq1k(q1)W(q2) - Dxq\ - D2ql)t\exp[iq2y] (5.187) xexp[iqi(x - sC{Ke{0) + Ke{qi))W{q2)t)}. The imaginary term is proportional to the velocity of the plane wave. Note that the real part of the linear growth rate is independent of the even part of the interaction kernel, while the velocity of the excitation is independent of the corresponding odd part. The instability criterion ReX > 0 can be written in the form: *<*>*<*> > (w+ if) • Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 110 where the notation for the ID bifurcation function A (5.176) is used. As mentioned, we consider only the case that diffusion in the y direction is much slower than that in the x direction, and put D<i = 0: WMKqi) > | j . (5.188) (Later I discuss qualitatively what happens at non-zero diffusion in the y direction.) The cases of "continuous" and "finite-range" interactions in the y direction are quali-tatively different and will be considered separately. The "continuous" case corresponds to the choice W = Wi, where the function W\ is given by (5.161) and has a positive Fourier transform (5.166). The main difference in the "finite-range" case (corresponding to the choice W = Wi-, see (5.161)) is that the function Wi has an oscillating Fourier transform (5.166). Recall that the functions Witi(qi) are even, so that the primary bifurcation always corresponds to iq^cj where qic is some critical wave vector. (i) T h e "continuous" case. In this case, the function W\{CJI) is positive and has a maximum equal to 1 at qi = 0 (see expression (5.166)). Then it follows from (5.188) that if instability occurs, it happens first at q2 = 0 (when the LHS of the instability criterion (5.188) is maximal). Then inequality (5.188) becomes and coincides with the ID instability criterion (5.178). This fact immediately results in the following predictions about the 2D "continuous" situation based on the ID results. Below, as in the ID case, I analyze conclusions from the instability criterion (5.188) at governing parameter values from different regions. Regions 1,2,10 (B < A,b> a{B/A)ll2): The 2D homogeneous steady state is stable because the bifurcation function is not positive, and the instability criterion (5.188) cannot be valid. Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 111 Regions 3,4,5,6 {B < A,b < a(B/A)^2 and B > A,b < a(B/A)^4): When diffusion is sufficiently slow, bifurcation caused by the mode qic ~ 0, q2c — 0 takes place. The growing pattern consists of long-wave harmonics in the x direction and is homogeneous in the y direction. Regions 7,8,9 (B > A,b > a(B/A)^4): Bifurcat ion from the homogeneous state occurs at q\c ^ 0,q2c — 0. The non-linear analysis in this case is completely analogous to the one in ID (Section 5.4.2). A pattern of "rolls" (stripes) emerges, with the axis of a stripe normal to the axis of a cell (Figure 24a). Though this pattern resembles rippling, the observed orientation of the stripes is different. The results are summarized in Table II. (ii) T h e "finite-range" case. The function W2(q2) n a s a n absolute maximum equal to 1 at q2 = 0 and an absolute negative minimum at q2c — 4.5. Negativity of the function W2(q2) at some values of the argument changes the bifurcation scenario radically. Again, the analysis is based on the properties of the bifurcation function described in Appendix D. Note: For the analysis in this section, we will need one more subdivision of the pa-rameter space with the dashed line in Figure 22. This line is defined by the trancendental equation: W2(q2c)\(0) = \(qlc), (5.189) where q2c ~ 4.5 and qic is given by (5.180). The solution of this equation cannot be obtained in an explicit form. However, it is not difficult to find the qualitative form of the line. We represent the equation of this line by b = WAA(B). Regions 1,10 (B < A,b > a(BfA)1^4): The bifurcation function X(qi) is negative, so at values of q2 for which W2(q2) is positive, there is no bifurcation. However, if the function ^2(92) becomes negative, the RHS of inequality (5.188) becomes positive. At Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 112 a sufficiently small value of the diffusion coefficient, a bifurcation can occur. Evidently, this happens first when the function W2(q2) reaches its absolute minimum, at the critical wave vector q2c — 4.5. The critical wave vector q\c should correspond to the absolute negative minimum of the function A(qi), which occurs at qic = 0. Thus, in the present case, the mode qic = 0,q2c — 4.5 breaks stability. This bifurcation corresponds to patterns of stripes parallel to the x axis, so that the axis of a cell is parallel to the stripe (Figure 24b). The non-linear analysis of such a pattern is very difficult. A multitude of neutrally stable non-periodic stationary patterns emerges after bifurcation. The pattern that actually develops depends upon the initial perturbation. This situation is as yet incompletely understood. Region 9 (and part of region 8 to the left of the dashed line) (B > A,b > WA,a{B)): The bifurcation function has a positive maximum at some finite qi, but it also has an absolute negative minimum at q\ = 0. The absolute maximum of the LHS of inequality (5.188) again corresponds to the values q\c = 0, q2c — 4.5 (see Appendix D). The character of the bifurcation and the non-linear analysis is the same as in the previous case. Thus, strong repulsion in the "finite-range" 2D case leads to some irregular behavior which is not fully understood. Region 7 (and part of region 8 to the right of the dashed line) (B > A, a(BfA)1/4 < b < WAta(B)): The maximum value of the bifurcation function is positive and occurs at q = qlc jL 0. In the region of model parameters to the right of the line b — WA,O.{B) defined in (5.189), the maximum of the LHS of the criterion (5.188) and, thus, the bifurcation, corresponds to the values q lc ^ 0 and q2c = 0 (see Appendix D). The patterns are growing "rolls" with bacteria normal to the axis of the roll. The non-linear analysis is similar to the ID case. Thus, in the 2D case of "finite-range" interaction, short-range attraction, long-range repulsion, a striped pattern occurs (Figure 24a). Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 113 Regions .{,5,6 (and part of the region 3 to the right of the dashed line) (B > A,b < a^B/A)1'4 and B < A,b < WA,O,{B)): The maximum value of the bifurcation function is positive and occurs at q = q\c = 0. The absolute maximum of the LHS of inequality (5.188) corresponds to q = q\c ~ 0 and q = qic = 0. This is the case of strong attraction, and long wave-length perturbations grow. Region 2 (and part of region 3 to the left of the dashed line) (B < A,WA,O,{B) < b < a(BI'A)1•' ): This case is the most interesting one. The bifurcation function is negative, so that bifurcation happens first at q — q2c = 4.5. The break of stability then corresponds to the negative minimum of the bifurcation function at q — q\c ^ 0 (qic is given by (5.189)). The leading unstable mode is characterized by the wave vector {(lie 7^ 0, qic = 4.5). Patterns of stripes, at some angle to the cell orientation other than 0 or 7r/2 are possible (Figure 24c). This pattern is precisely the case of the rippling phenomenon. The conclusion is that this particular phenomenon takes place when the following conditions are met: (i) There is "finite-range" interaction. (ii) There is an appropriate balance between the short-range repulsion and long-range attraction. The results of linear stability analysis are summarized in Table 2. Recall that due to the symmetric properties of the interaction kernels, the critical wave vectors in the case of rippling are (±qi c 7^ 0, ±q2C 7^ 0). Then another possibility is a rhombic (spots) pattern (Figure 24d) which also resembles rippling. To distinguish between these two possibilities, non-linear analysis will be carried out in the next subsection. If D2 7^ 0, the results concerning the case of "continuous" interaction does not change qualitatively. Also, in the case of "finite-range" interaction, if the values of interaction Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 114 parameters are in regions 2-7, the anaysis presented here gives qualitatively similar re-sults. However, if the model parameters are located in regions 1,8-10, an "orientational" transition occurs as diffusion in the y direction grows. The rolls shown in Figure 24b "turn" and become stripes similar to those in Figure 24c or rombs as in Figure 24d. This fact implies the possibility of a diffusion-driven bifurcation to the rippling structure. The biological meaning of this result is as yet unclear. Before concluding the linear stability analysis, consider the imaginary term in the dispersion relation (5.186). This term indicates that growing patterns are travelling waves. (It is easy to see that standing waves are excluded.) The absolute value of this term defines the velocity of the "ripples". Since the coordinate system moves with the velocity of an individual bacterium, the non-zero imaginary term in (5.186) implies that the velocity of the ripples may differ from the individual velocity as well as from the velocity of the swarm as a whole. Non-linear analysis for s ^ 0 is very cumbersome. At present, we describe the non-linear analysis without this term. 5.5.2 2 D Non-l inear analysis I perform here non-linear analysis of the 2D Model (5.162) in the limiting case s = 0. The general case would be qualitatively similar, but much more cumbersome. The approach here is similar to the ID case discussed in section 5.4.2. We expand the density distribution as follows: C = I + Gi cos (qlcx + q2cy) + G2 cos (qlcx -q2cy) + ^ d(<lu 92) cos (q^ + q2y), (9l,92)^(gic.±92c) (5.190) where G\ and G2 are amplitudes of the unstable modes (characterized by wave vectors (qic-, ±Q2c) respectively. Denote qc = (q lc, q2c). Then g(qi-,q2) are amplitudes of the stable modes. Among the stable modes, dominating ones are those characterised by the wave Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 115 vectors q = (2qic, 0) and q = (2qic, 2q2c). The asymptotic equations for the corresponding amplitudes are: 9(q) = A(q) g{q) - qicR(qc)Gl, g(q) = \(q)g(q) - qlcR(qc)GxG2, \(q) = \Dxq\c + AD2q22c + 2qlcR(q), X{q)=AD1q2le + 2qlcR(q)i R{qc) = K{qu)W2{q2c), R(q) = K\2qlc)W2(2q2c), R{q) = K(2qlc). These equations were derived by plugging (5.190) into (5.162) and keeping leading terms corresponding to the wave vectors q and q. Assuming a stationary situation, mode am-plitudes in these equations can be expressed as quadratic polynomials of the amplitudes of the leading modes. The latter, on the other hand, depend on the amplitudes of the dominating stable modes as follows: G1 = 6\G1-pg(q)G1-ug(q)G2, where 1 1 P = 2<?ic {Hie) + R{qj) , u = -qu (k{qc) + R{q)) , SX = Dxq\c + D2q\c + qlcR(qc), (and a similar equation for G2. The present approach is valid close to the bifurcation, when 5\ is small and positive.) Solving the equations for the stable dominant ampli-tudes and plugging the resulting algebraic expressions into the equations for the unstable amplitudes, we finally get the following closed system of generalized Ginsburg-Landau ODEs: Gi = 8\Gi - pG3, - uGlGi, 1 2 (5.191) G2 = SXG2-pGl-uGlG2, Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 116 w here P = jrk-Hqc) {Hqc) + h{qj), (5.192) .2 u = -^R(qc) (k{qe) + % ) ) • (5.193) In the region of interest, and with the present choice of the interaction kernels , the Landau constants p and u are positive. If for some other choice of kernels these constants are negative, the present analysis is invalid, and the exact behavior is still an open question. The system (5.191) allows three qualitatively different steady states: 1) Gi = G2 = 0, 2) Gi = 0, G2 = J—, 3) GX = G2 = \ / - ^ - . (5.194) V p V P + U The first solution corresponds to the homogeneous steady state and is unstable. The second state describes a striped pattern, and the third describes a spotted pattern. The leading mode amplitudes in the vicinity of the striped steady state are governed by the linearized system of ODE's: Gi = 8\ (1 — -] G\, u> p — stripes stable, 2) I / V v) (5.195) Gi = — 25\Gi, u < p — stripes unstable. The striped pattern is locally stable if the condition u > p is met (then the eigenvalues of the stability matrix of the system (5.195) are then strictly negative). If the opposite inequality occurs, this structure is unstable. In the vicinity of the spotted steady state the leading amplitudes obey the equations: G\ = {pG\ -f uG2) u > p — spots unstable, p + u , 3) P28x (5.196) (?2 = {uG\ + pG2) u < p — spots stable. I p + u The spotted pattern is locally stable if the striped structure is unstable (if the condition u < p is met then the eigenvalues of the stability matrix of the system (5.196) are strictly Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 117 negative). Conversely, stripes are stable if spots are not. Thus, in order to understand which pattern is realized, numerical estimates of the values of the Landau constants are needed. As an example, the following values of the parameters of the 2D model were used: A = 4, B = l, a = 2, 6 = 1 , C = 1, D2 = 0. Using the formulae derived in this chapter and the program Mathematica, the values of the critical wave vectors were found to be q\c = v2,q%c — 4.5, the critical diffusion coefficient was D\ ~ 0.15, and the Landau constants were: u ~ 0.07 > p ~ 0.035. The conclusion is that the striped pattern is stable at these parameter values. The possibility of a spotted (rhombic) pattern at some other choice of interaction kernel is not ruled out. 5.6 Discuss ion The main result obtained in this chapter is the existence of a specific pattern of parallel stripes. The pattern is such that the angle between the cell axes within a stripe and the stripe axis is different from 90°. The conclusion from the linear stability analysis and weakly non-linear analysis is that the appearance of this pattern is very sensitive to the character of the interactions between cells, and requires significant long-range re-pulsion between cells. In the case of chemotactic interaction, both chemoattractant and chemorepellent should be present, and a certain very delicate balance has to be main-tained between the two substances. This pattern resembles the experimentally observed periodic trains of ripples in a Myxobacteria swarm. Chapter 5. Partial integro-differential model for aggregation in a bacterial swarm 118 The existence of other types of patterns, which are not observed in Myxobacteria are also predicted. Their description may be useful for the analysis of other types of aggre-gation phenomena. At a certain balance between cell-to-cell repulsion and attraction, the striped pattern that appears has stripes normal to the direction of motion of the swarm. Significant attraction between bacteria would lead to the creation of a compact, dense, unstructured swarm. A peculiar phenomenon may arise, when a periodic pattern turns into a dense aperiodic swarm as swarm density grows. In theory, the appearance of rhombic (spot) patterns is also possible. Emerging patterns depend in a sensitive and non-monotonic way on diffusion in a direction normal to the cell's orientation. Some problems and questions are still open. I performed non-linear analysis in the case that the even convolution term is absent. This term would be responsible for the density-dependent average net drift of the cells. Some preliminary investigation shows that the growth of this neglected term suppresses periodic patterns. The explanation and biological meaning of this behavior is still open. I analysed the development of instabilities of the homogeneous pattern. Similar analysis is desirable for stability break from the travelling wave pattern. What is the shape of the swarm edge then? Which patterns emerge when the weakly non-linear analysis fails? There may be non-stationary solutions which are chaotic in some sense (D. Wollkind, private communication). Possibly, the mathematical problem may be ill-posed. Chapter 6 Discuss ion 6.1 Resu l t s In Chapter 2 we explored three models for self-alignment. Though in each model the specific mechanisms of alignment were different, the models are united by the following common feature. Cells or molecules turn randomly and independently, while the interac-tions lead to one or other type of alignment. In other words, there is a superposition of linear terms, representing dispersal and random motion and nonlinear terms, represent-ing interaction. As either the density of the population grows, or the rotational diffusion slows down, the interaction terms prevail, causing a bifurcation. It was shown that close to the bifurcation, the homogeneous angular distribution was slightly perturbed by one leading mode. The form of this mode depends on certain symmetry properties of the integral kernels, rather than on their detailed form. Biologically, the results are simple and intuitively clear. If "head-to-head" alignment is preferable (in a probabilistic sense), then a unipolar bundle starts to form. If "head-to-head" and "head-to-tail" alignment are equally likely, then a bipolar bundle develops (in the simplest case, there is no difference between the "head" and the "tail"). Finally, if interaction is most favorable at a relative angle of 90°, then orthogonal networks evolve. These results obtained in Chapter 2 for values of the model parameters close to the bifurcation are complemented by the results of Chapter 3 for values far from the 119 Chapter 6. Discussion 120 bifurcation. We have shown that if the rotational diffusion becomes very slow, or the density very high, then mild inhomogeneities in the angular distribution become more prominent, and in the zero diffusion limit evolve into narrow peaks. The stability of these peaks is not proved rigorously, but is well supported by the analytical treatment of certain class of perturbations and by numerical experiments. The width of these peaks was estimated and found to be small (of order e a , where e is a small renormalized rotational diffusion coefficient and values of a were found to be equal to 1/2 and 1/4 for Models II and I, respectively). The biological meaning of these results is, again, easy to understand. When the tendency for disorder decreases or the concentration grows (effectively promoting alignment), then the steady state approaches the phase of almost complete alignment. In Chapter 4, linear stability analysis demonstrated the existence of two kinds of bifurcations distinct from the purely angular one. The first one describes the formation of dense angularly disordered patches. This happens when random rotation of the objects is relatively fast (tipping the balance away from angular order), while an attraction in space dominates over the spatial dispersal. Another bifurcation, the most complex one, signifies either an alignment correlated in space or patches formed by aligned individuals. The character of the macroscopic pattern evolved depends on the microscopic details of the pair-wise interaction between individuals as a function of their mutual distance and orientation. The important feature of the spatial bifurcation is the smallness of wave vectors responsible for the bifurcation, implying a long length scale of the spatial pattern near criticality. This means that the length scale of the pattern may be comparable to the size of the spatial domain and the effect of the boundaries could be significant. One of the implications is that irregular "droplets" containing large numbers of objects start to Chapter 6. Discussion 121 grow. This is similar to patterns formed in fibroblast cultures. The main result of Chapter 5 is that certain simple assumptions about Myxobacteria interactions are sufficient to account for the peculiar "rippling" pattern of stripes in the density profile of the swarm. The dependence of the interactions on cell-cell distances is important. The interaction may be either "direct", through some sort of signalling, or "indirect", through chemoattractants and chemorepellents. In the latter case, the presence of both types of chemicals is necessary. In both cases, the interaction should be repulsive at short distances and should decrease in a non-smooth fashion in the direction normal to the gliding motion in order for the rippling phenomenon to occur. 6.1.1 Significance of the results The proposed type of modelling has several advantages. First, the consideration of joint spatio-angular distributions, elucidates the interconnection between angular and trans-lational modes of movement. Second, non-local effects are treated in a more complete way than before without reduction to PDEs. This analysis leads to an explanation of short-range long-lasting inhomogeneities. Finally, the mathematical description of pat-terns was extended to the more realistic three-dimensional cases. Apart from the purely mathematical significance of the research done for the theory of non-linear partial integro-differential equations, applications to biology and biophysics hold some promise. 6.2 O p e n quest ions One of the questions which is left open is the temporal behavior of the biological systems in question. How fast do the predicted spatio-angular patterns evolve from an arbitrary initial state? What are the characteristic relaxation times of the system? In other words, how fast do small perturbations of the stable pattern die out? A quantitative answer to Chapter 6. Discussion 122 the first question can be given based on numerical simulations. The second problem is treatable analytically in the limiting case of slow diffusion. Qualitatively, in this case, objects wander in an "external field of forces" created by the stable non-homogeneous density. The fact that the deviation from this stable pattern is small means that interaction between the perturbed objects can be neglected. Mathe-matically this situation is modelled by a linear Fokker-Planck equation. Initial analysis leads to the result that the relaxation time is exponentially long in this case. The im-portant biological implication of this fact is that, in reality, relaxation times may be comparable to the existence time of the system in a stable environment. Thus, a vast variety of quasi-stable long-lived patterns not predicted by the model may occur. In the limiting case of almost complete alignment, some other angular distributions (either smooth, or peak-like, but with other configurations of peaks) may evolve. The problem of an analytical approach to the corresponding stability analysis would lead to complicated integro-differential equations. It would be also of great interest to undertake a nonlinear stability analysis of the spatio-angular distributions in order to find out the type of order actually evolving. Another possible approach to this problem is to try to derive some simpler toy models allowing explicit solutions (such as, for instance, the aggregation model of Grinrod (1991)). The models in Chapters 2 - 4 capture some of the important features of aggregation and alignment phenomena, but do not describe the details of these phenomena. It seems that more specific equations modelling, for example, fibroblast cultures in detail would be resistant to analytical treatment. A reason for this pessimism is that , most probably, linear operators of a realistic model would not share the same complete and orthogonal system of eigenfunctions. In that case even a linear stability analysis would become too complex to be informative. Nevertheless, this more biologically oriented modelling should Chapter 6. Discussion 123 be undertaken. The emphasis then would be on numerical methods. The most interesting problem stemming from Chapter 5 is to design some biological experiments in order to check predictions of the model. Arranging different densities of the swarm under otherwise equal conditions would be one such experiment. Preliminary discussions with D. R. Zusman, though, indicate, that such experiments are not very easy to arrange. On the mathematical side, one could speculate that for some choice of the interaction kernels families of explicit solutions in the form of trains of travelling waves can be obtained. It would be worth finding these kernels and solutions, and investigating their stability. Also, a more realistic model of the bacterial swarm, taking into account turning modes of motion should be considered. This would be similar to the one introduced in Pfistner and Alt (1989). A numerical treatment of the nonlinear hyperbolic equations of Chapter 5 and even simpler cellular automata experiments promise to provide valuable insight. Finally, the Lagrangian approach to modelling bacterial swarms should be successful in limiting situations of very high swarm density or very slow spatial diffusion. Then, as in Model II in Chapter 3, one can demonstrate that the integro-differential equa-tions of Chapter 5 can be written in a gradient form. A qualitative investigation of the corresponding Lyapunov function allows one to draw conclusions about limiting spatial patterns under assumptions about parameters of the model. In particular, it can be shown that if spatial attraction prevails, the swarm "collapses" and a single very dense patch of bacteria forms. This indicates the limitations of our model: anti-crowding mech-anisms should prevent such "collapse". With such mechanisms included explicitly one would expect the development of a single unstructured bacterial patch. Chapter 6. Discussion 124 6.3 S o m e o t h e r app l i ca t i ons In this section I mention some related phenomena which were beyond the scope of this Thesis. Alignment phenomena occur often in nonbiological systems. In liquid crystals, rod-like molecules form various associations, including nematics, smectics, and cholester-ics (de Gennes, 1974; Chandrasekhar et al., 1970; Greco and Marrucci, 1992; Luckhurst and Gray, 1979; Sheng, 1973). The type most closely related to the kind of order investi-gated in Chapters 2 and 3 are the nematic liquid crystals in which the molecules all have a common axis of alignment, but no long-range spatial order such as a regular crystal structure. In Chapter 4, where both spatial and angular order are investigated, we find analogies with smectic-like order, where both long-ranged spatial order and alignment occur. We briefly comment on other applications of integro-partial differential equations in biological modelling. A review of integral equations in biology is given in Levin and Segel (1985). One body of theory that has always relied heavily on the formulation and analysis of integral equations is the theory of activity and evolution of networks of biological neurons. An early example of this type of modelling is to be found in papers by Swindale (1980, 1991) who describes development of patterns of ocular dominance or orientational response in cells of the visual cortex of mammals. Similarly, models that concern neural interactions and thus contain convolutions appear in Ermentrout and Cowan (1979), Ermentrout et al. (1986). In such models, the activation of a neuron at site x is, in general, dependent on inputs and stimuli arriving from neurons at remote sites, whose axons impinge on and synapse with the given local nerve cell. The type of effect of one neuron on another may typically depend on the mutual distance of the pair, and it is customary to assume that nearby neurons activate one another, whereas more distant ones have mutually inhibitory influences. For this reason, equations containing Chapter 6. Discussion 125 integral terms (which sum up the contributions of all neurons) with convolutions (that describe how the signal depends on the mutual distance) are convenient representations of these phenomena. Such equations do not come from balance arguments, and indeed the quantity of interest (which may be the intensity of activity) is not generally a conserved quantity. Many organisms do not move in a smooth continuous way, but rather execute a series of small bursts or discrete jumps that can have a variety of possible sizes. This type of motion has been modelled in the traditional equations of mass balance by inclusion of integral terms with kernels that represent the jump size distribution. Examples of this approach include Alt (1988), Othmer et al. (1988). An application of integral equations to density of branches in a network of filaments that interact by cross-linking was inves-tigated by Edelstein-Keshet and Ermentrout (1989). Models in epidemiology deal with temporal and spatial distributions of infectuous and susceptible individuals where long-range effects are important (Kot, 1992). Finally, evolution of patterns in aspect (diver-sity of colors among prey species; assortative mating with distributed density-dependent rates; dominance in bumble bee colonies) can be described in terms of dynamics in some abstract aspect space with a long ranged interaction Levin and Segel (1985). Integral equations have also appeared recently in the literature on coupled oscillators. An evolution equation for the density of oscillators at phase 9 and frequency to, analogous to a model by Kuramoto (1975) was derived and analyzed by Strogatz and Mirollo (1991). The oscillators interact with one another with intensity that depends on their relative phases, and this causes changes in the frequency and phase. Trains of travelling wave patterns are well-known and widely studied in hydrodynam-ics. Integro-PDE's were also used for this goal (Benjamin, 1974). Stripes and rhombic convective patterns were discussed in Malomed and Tribelskii (1987). Travelling wave Chapter 6. Discussion 126 train solutions caused by long-range effects have been discussed in this Thesis. Such so-lutions without non-local interactions can take place for reaction-diffusion systems with limit cycle kinetics (Howard and Kopell, 1977). These waves appear under certain con-ditions in excitable media with multiple steady states and diffusion. Such waves have various features, and may include slowly varying waves, shock structures etc. These may successfully describe a range of dispersal patterns found in nature. A periodic peak-like ID pattern in a chemotactic system was obtained numerically in Myersough and Murray (1992). Finally, it is worth mentioning that similar phenomena occur in some applications outside the natural sciences. In traffic flow, interaction between cars causes instability of homogeneous traffic. The equations derived in the continuous approximation are quite similar to the ones describing a propagating swarming and dispersing biological popu-lation (Kerner and Konhauser, 1993). The existence of travelling waves was shown in Kerner and Konhauser (1993). Long ranged effects, which are probably important, were not yet taken into account in any of these traffic flow models. 6.4 S o m e final notes The usual approach to the theoretical description of aggregation and alignment patterns in physics is to introduce a free energy functional and to treat the pattern as the function minimizing this functional. As we have seen, an analogous technique works in some limiting cases of our models. There is no reason, however, that in a highly non-equilibrium biological phenomenon, the same approach would succeed. On the other hand, our phenomenological approach has the flaw of being too de-tailed. Theoretical physics experience teaches that if the system is highly complex, then it is best to work with the so-called order parameters, rather than with detailed spatio-angular distributions. Mathematically that would lead to a significant reduction of the Chapter 6. Discussion 127 dimensionality of the problem, and would lead to the so-called Ginzburg-Landau equa-tions. There is vast experience from physics available to treat such equations, and this direction remains open for future exploration. One of the peculiar conclusions from our investigation is that the character of a pat tern evolving after the bifurcation is not sensitive to the detailed form of the kernels in the angular case (Chapters 2 - 4 ) , but it is very sensitive in the spatial aggregation case (Chapter 5). We can give a speculative explanation of this phenomenon. Angular space is compact (and the corresponding kernels expand into Fourier series), while real space is not compact (and the kernels expand into Fourier integrals). In a way, this means that many more modes compete in the spatial case, leading to a greater variety of patterns. We have left untouched the question of well-posedness of the kind of integro-differential equations that have appeared in this Thesis. This remains a grave problem. It is known that the description of aggregation with "anti-diffusion" is ill-posed, but the introduction of non-local effects in the form of integral terms makes the problem well-posed. This is one more open area of research. Bibliography ABRAMOWITZ, M. AND STEGUN, L A . 1970. Handbook of Mathematical functions with formulas, graphs, and mathematical tables. National Bureau of Standards, Ap-plied Mathematical Series. 55. ALT, W . 1985. Degenerate diffusion equations with drift functionals modelling ag-gregation. Nonlin. Analysis 9, 811-836. ALT, W . 1988. Modelling of motility in biological systems. ICIAM 1987 Proceedings, SIAM Philadelphia, 15-30. ARNOLD, V . I. 1978. Mathematical methods of classical mechanics. Springer, NY. BEHMLANDER, R . M . AND DWORKIN, M. 1994. Biochemical and structural anal-yses of the extracellular matrix fibrils of Myxococcus xanthus. J. Bacteriol. 176 , 6295-6303. BENJAMIN, T . B . 1974. Lectures on non-linear wave motion. Non-linear wave mo-tion. Ed.: Newell, A. C , ed., AMS, Providence, 3-47. BEN-YAKOV, E. , SHOHET O.AND T E N E N B A U M A. 1994. Generic modelling of cooperative growth patterns in bacterial colonies. Nature 368 , 46-49. BLATT, Y U . AND ElSENBACH, M. 1995. Tar-dependent and Tar-independent pat-tern formation by Salmonella typhimurium. manuscript. BRAY, D . 1992. Cell movements. Garland Publ., NY. BUDRENE, E . O . AND BERG, H.S. 1991. Complex patterns formed by motile cells of E.coli. Nature 349 , 630-633. BUSSE, F . H . 1987. Patterns of bifurcation from spherically symmetric cases. The Physics of Structure Formation. Ed.: W. Guttinger, G. Dangelmayr, Springer Ver-lag, NY, 88-96. CHANDRASEKHAR, S., SHASHIHAD, R.AND TARA, N. 1970. Theory of melting of molecular crystals: the liquid crystalline phase. Molec. Cryst. Liq. Cryst. 10, 337-358. 128 Bibliography 129 [13] CHANG, K. W . AND HOWES, F . A. 1984. Nonlinear Singular Perturbation Phe-nomena: Theory and Application. Springer, NY. [14] ClVELEKOGLU, G. AND EDELSTEIN-KESHET, L. 1994. Models for the formation of Actin Structures. Bull Math Biol 56, 587-616. [15] DE GENNES, P . G. 1974. The Physics of Liquid Crystals. Claredon Press, Oxford. [16] DWORKIN, M. AND KAISER, D . 1985. Cell interactions in myxobacterial growth and development. Science 230 , 18-24. [17] EDELSTEIN-KESHET, L. AND ERMENTROUT, G. B. 1989. Models for branching networks in two dimensions. SIAM J Appl Math 49 (4) , 1136-1157. [18] E D E L S T E I N - K E S H E T , L. AND E R M E N T R O U T . G .B . 1990. Models for contact me-diated pattern formation. J. Math. Biol. 29, 33-58. [19] ELSDALE, T . 1972. Pattern formation in fibroblast cultures, an inherently precise morphogenetic process. Towards a Theoretical Biology Ed.: C. H. Waddington, Edinburgh University Press, Edinburg, 22-36. [20] ELSDALE, T . AND W A S O F F , F . 1976. Fibroblast cultures and dermatoglyphics: the topology of two planar patterns. Wilhelm Roux's Archives 180, 121-147. [21] ERMENTROUT, G. B., CAMPBELL, J. AND OSTER, G. 1986. A model for shell patterns based on neural activity. The Veliger 28(4) , 369-388. [22] ERMENTROUT, G. B. AND COWAN, J. 1979. A mathematical theory of visual hallucination patterns. Biol. Cyber. 34, 137-150. [23] E R M E N T R O U T , G. B . AND E D E L S T E I N - K E S H E T , L. 1991. Cellular Automata ap-proaches to biological modelling. J. theor. Biol. 160, 97-133. [24] ERMENTROUT, G. B . AND M C L E O D , M. C 1993. Existence and uniqueness of travelling wave for a newral network. Proc.Roy.Soc.Edinburgh 123A , 461-478. [25] FRIEDRICH, R. AND H A K E N , H. 1989. A short course on synergetics. Nonlinear phenomena in complex systems. Ed.: A.N. Proto, Elsevier Publ., North Holland, 103-150. [26] GANTMAKHER, F . R. 1966. Theory of Matrices. Springer, NY. [27] GARDINER, C. W . 1985. Handbook of Stochastic Methods. Springer, NY. Bibliography 130 [28] GRECO, F . AND MARRUCCI, G. 1992. Rodlike molecular dynamics, the tumbling regime. Mol. Crystals Liq. Crystals 212, 125-137 [29] GRINROD, P . 1991. Patterns and Waves. Claredon Press, Oxford. [30] GRINROD, P . 1988. Models of individual aggregation in single and multi-species communities. J. Math. Biol. 26, 651-660. [31] GRINROD, P . , SlNHA, S. AND MURRAY, J . D. 1989. Steady-state patterns in a cell-chemotaxis model. IMA J. Math Appl in Med. Biol. 6, 69-79. [32] GROSS, M. C. AND HOHENBERG, P . , C . 1993. Pat tern formation outside of equi-librium. Rev. Mod. Phys. 65 , 851-112. [33] GRUNBAUM, D. 1994. Swarming behaviour as an aide to Chemotaxis. 3D Animal Aggregations. Ed.: J. Parrish, W. Hammner, Cambridge Univ. Press, 1-35. [34] GRUNBAUM, D. AND OKUBO, A. 1994. Modelling social animal aggregation. Fron-tiers in Mathematical Biology. Ed.: S. Levin, Springer, NY, 296-325. [35] HAKEN, H. 1977. Synergetics. Springer Verlag, New York. [36] H O L M E S , E . E . , L E W I S M. A., B A N K S J . E . AND V E I T R .R . 1994. PDE in ecology: spatial interactions and population dynamics. Ecology 75(1) , 17-29. [37] HOWARD, L . N . AND K O P P E L , N . 1977. Slowly varying waves and shock structures in reaction-diffusion equations. Studies in Appl. Math. 56, 95-145. [38] HUNDING, A. 1982. Spontaneous biological pattern formation in 3D sphere. Prepat-terns in mitosis and cytokinesis. Evolution of Order and Chaos. Ed.: H.Haken, Springer, Berlin, 100-128. [39] HUTH, A. AND WiSSEL, C. 1989. The movement of fish schools: a simulation model. Biological Motion. Ed.: W. Alt, G. Hoffman, Springer, NY, 577-590. [40] JAGER, E . AND SEGEL, L. 1992. On the distribution of dominance in population of social organisms. SIAM J. App. Math. 52-5, 1442-1468. [41] KAISER, D . AND CROSBY, C. 1983. Cell movement and its coordination in swarms of Myxococcus xanthus. Cell Motil. 3 , 227-245. [42] KELLER, E . F . AND SEGEL, L. A. 1971. Travelling bands of chemotactic bacteria: a theoretical analysis. J. Theor. Biol. 30, 235-248. Bibliography 131 KERNER, B . S. AND KONHAUSER, P . 1993. Cluster effect in initially homogeneous traffic flow. Phys. Rev. E48, R2335-2338. KOT, M. 1992. Discrete-time travelling waves: ecological examples. J. Math. Biol. 30, 413-436. KRAUT, E. A. 1979. Fundamentals of Mathematical Physics. Kreiger Publ, NY. KURAMOTO, Y. 1975. Self-entrainment of a population of coupled nonlinear oscil-lators. International Symposium on Mathematical Problems in Theoretical Physics. Ed.: H. Araki, Springer, NY, 420-422. LANDAU, L. D . AND LiFSHlTZ, E . M. 1956. Hydrodynamics. Pergamon Press, London. LANDAU, L. D . AND LiFSHlTZ, E . M. 1976. Mechanics Pergamon Press, Oxford. LAUFFENBURGER, D . A. 1984. Chemotaxis and cell aggregation. Modelling of pat-terns in space and time. Ed.: W. Jager, J. D. Murray, Springer, Berlin, 197-213. LEVIN, S. A. AND SEGEL, L. A. 1985. Pat tern generation in space and aspect. SIAM Review 27, 45-67. LEWIS, M. A. 1994. Spatial coupling of plant and herbivore dynamics. J. Math. Biol. 45 , 277-312. LiFSHlTZ, E . M. AND PiTAEVSKlI, L. P . 1981. Physical Kinetics. Pergamon Press, Oxford. LOGAN, J . D . 1994. An introduction to nonlinear partial differential equations. J. Wiley and Sons, NY. LUCKHURST, G. R. AND GRAY, G. W . 1979. The Molecular Physics of Liquid Crystals. Acad. Press, NY. MACROBERT, T . M. 1967. Spherical Harmonics. Pergamon Press, Oxford. MARCUS, M. AND MING, H. 1964. A Survey of Matrix Theory and Matrix Inequal-ities. Align and Bacon, Boston. [57] MALOMED, B . A. AND TRIBELSKII, M. I. 1987. Stability of stationary periodic structures for weakly supercritical convection and in related problems. Sov. Phys. JETP 65 , 305-310. Bibliography 132 MATSUYAMA, T . AND MATSUSHITA, M. 1993. Fractal morphogenesis by a bacterial cell population. Critical Reviews in Microbiol. 19, 117-135. MOGILNER, A. AND EDELSTEIN-KESHET, L. 1995a. Selecting a common direc-tion I: How orientational order can arise from simple contact responses between interacting cells. J. Math. Biol. 3 3 , 619-660. MOGILNER, A. , EDELSTEIN-KESHET, L. AND ERMENTROUT, G. B . 1995b. Se-lecting a common direction ILPeak-like solutions representing total alignment of cell clusters. Submitted for publication. MOGILNER, A. AND EDELSTEIN-KESHET, L. 1995c. Spatio-angular order in popu-lations of self-aligning objects: formation of oriented patches. To appear in PhysicaD. MORSE, P . M. AND FESHBACH, H. 1953. Methods of theoretical physics. McGraw-Hill, NY. MURRAY, J . D . 1989. Mathematical Biology. Springer, Berlin. MYERSCOUGH, M.R . AND MURRAY, J . D . 1992. Analysis of propagating pat tern in a chemotactic system. Bull.Math.Biol. 54, 77-94. NAGAI, T . AND MlMURA, M. 1983. Asymptotic behavior for a non-linear degener-ate diffusion equation in population dynamics. SIAM J. Appl. Math. 43 , 449-464. O K U B O , A. 1990. Diffusion and Ecological Models. Springer, NY. O'MALLEY, R. E . J R 1974. Introduction to Singular Perturbations. Academic Press, NY. ONSAGER, L. 1949. The effects of shape on the interaction of colloidal particles. Ann. NY Acad. Sci. 5 1 , 627-659. OTHMER, H. G., DUNBAR, S.R. AND ALT, W. 1988. Models of dispersal in biological systems. J. Math. Biol. 26, 263-298. PENROSE, O. 1978. Kinetics of phase transitions. Stochastic Processes in Nonlinear Systems. Ed.: L. Garrido, Springer, NY, 211-222. PFISTNER, B . 1989. A one dimensional model for the swarming behavior of myxobacteria. Biological Motion. Ed.: W. Alt, G. Hoffmann, Springer, NY, 556-563. Bibliography 133 PFISTNER, B . AND ALT, W . 1989. A two-dimensional random walk model for swarming behaviour. Biological Motion. Ed.: W. Alt, G. Hoffmann, Springer, NY, 564-565. PRIESTLY, E. B., WOJTOWICZ, P. J. AND CHENG, P. 1975. Introduction to Liquid Crystals. Plenum Press, NY. REIHENBACH, H. AND DWORKIN, M. 1981. Introduction to the gliding bacteria. Procaryotes. Ed.: M. P. Starr, H. Stolp, Springer, Berlin, 113-192. REICHENBACH, H. 1984. Myxobacteria: a most peculiar group of social procaryotes. Myxobacteria: development and cell interaction. Ed.: E. Rosenberg, Springer, NY, 2-54. REICHENBACH, H. 1986. The Myxobacteria: common organisms with uncommon behavior. Microbiol.Sci. 3 , 268-274. ROSENFELD, L. 1948. Nuclear forces. North Holland Publ. Co., Amsterdam. SAN MIGUEL, M. AND SAGUES, F . 1990. Transient pattern Dynamics. Patterns, defects and material instabilities. Ed.: D. Walgraef, N. M. Ghoniem, Kluwer Publ., Dordrecht, 13-26. SEGEL, L. A. 1984. Modelling Dynamic Phenomena in Molecular and Cellular Biology. Cambridge University Press, Cambridge. SHAPIRO, J . A. AND TRUBATCH, D . 1991. Sequental events in bacterial colony morphogenesis. Physica D49 , 214-223. SHENG, P . 1973. Effects of bundling in a lattice gas model of liquid crystals. J. Chem. Phys. 59, 1942-1952. S H I , W . AND Z U S M A N , D . R. 1994. Sensory adaptation during negative chemotaxis in Myxococcus xanthus. J. Bacteriol. 176, 1517-1520. SHIMKETS, L. J . AND KAISER, D. 1982. Induction of coordinated movement of Myxococcus xanthus cells. J .Bacteriol. 152, 451-461. SHIMKETS, L . J . 1990. Social and developmental biology of Myxobacteria. Micro-biol.Reviews 54, 473-501. SMITH, D. R. 1987. Singular Perturbation Theory, an Introduction with Applica-tions. Springer, NY. Bibliography 134 [85 [86 [87. [88 [89 [90 [91 STROGATZ, S. H. AND MlROLLO, R. E . 1991. Stability of incoherence in a popu-lation of coupled oscillators. J. Stat. Phys. 63, 613-635. SwiNDALE, N. V. 1980. A model for the formation of ocular dominance stripes. Proc. R. Soc. Lond. B 208, 243-264. SwiNDALE, N. V. 1991. A model for the coordinated development of columnar systems in primate striate cortex. Biol. Cybern. 66, 217-230. TABONY, J. AND JOB, D . 1990. Spatial structures in microtubular solutions re-quiring a sustained energy source. Nature 346, 448-451. TURCHIN, P . 1986. Models for aggregating populations. Mathematical Ecology. Ed.: T. G. Hallam et al., World Scientific, Singapore, 101-127. ZWANZIG, R. 1963. First-order phase transitions in a gas of long thin rods. J. Chem. Phys. 39, 1714 - 1721. ZWILLINGER, D. 1990. Handbook of differential equations. Acad.Press, NY. A p p e n d i x A Bifurcations in an angular space A . l Operators and eigenfunctions in 3 D The Laplacian operator in surface spherical coordinates has the form: i f> fin i fan AC(9A,t) = _ -^ ( 8 in^) + -^s~- (A.197) s i n <p a<f> o<p s i n <p ad1 The convolutions expressed in terms of the angles <f> and 9 are K*C= / K{$,6,cj)',0')C{(t>',e')sm(l>'de'dcj)'. (A.198) Jo JO The Laplacian operator in surface spherical geometry has the eigenfunctions Yn(<f>,6) of degree n with AYn = -n{n + \)Yn. (A.199) Yn can be expressed in the form: n Yn(<f>,0) = A0P°(cos<£) + ^ ( ^ m C O s m 0 + B m s i n m ^ ) Pnm(cos<^), (A.200) m = l where A m ' s and B m ' s are arbitrary constants, P° are Legendre polynomials of degree n, P™ are associated Legendre functions of degree n and order m. (Macrobert, 1967; Kraut, 1979.) The first few Legendre polynomials for the case m = 0, P° , (usually written simply 135 Appendix A. Bifurcations in an angular space 136 Pn with the superscript omitted.) written in terms of the argument cos (f> are: < Po(cos^) = 1, Pi (cos 4>) = COS(j), P2(cos0) = §(cos</>)2- | . (A.201) The inner products of the SSH with the Legendre polynomials are also SSH as follows: Pn(n) y , (<M') d (cos^ ) (# ' = — — yn(<£,0) <Jn,,, (A.202) -l In + i 1, n = /; 0, ra^Z. (See Macrobert (1967) or any classical text on spherical harmonics.) This means that where where J K(Sl - SI') Yn{<j,',6') dtt = Kn Yn(<f>,9), Kn = K(n) = 2TT £ K{rj) P°(r,) dV. (A.203) (A.204) This establishes that the SSH are eigenfunctions of both the Laplacian and the operator K* in 3D. A.2 Stabil i ty of the homogeneous s teady s tate in 3 D In 3D, stability analysis of equations (2.20) to perturbations of the form (2.29) leads to a Jacobian matrix, J = —n —8/ (A.205) Appendix A. Bifurcations in an angular space 137 r2w Km = A 3cos[^ 7] , cos7>3 /4 ; (A.206) where the constants stand for the following 8 = /in(n + 1) + / ^ ( l + £ n ) + PP. Thus, we find that in 3D det J = ^£^[^-n(n + 1) - (^)2Kn(l - Kn)}. (A.207) p 7 7 Thus, in 3D, instability occurs whenever An(n + 1) < Kn(l - Kn), A = -{^~j)2. (A.208) 7 piw Example kernels The exact functional forms of the kernels used to produce Figure 3 were: ( Ai(cos7 — 3/4)-1/2, c o s 7 > 3 / 4 ; Ki = { 0, cos 7 < 3/4. A2, co s7>3 /4 ; Kn = { (A.209) 0, c o s 7 < 3 / 4 . 0, cos 7 < 3/4. High harmonics are stable Since A in equation (A.208) is positive, instability is most likely for low values of the integer n, e.g., n = 1 and high harmonics cannot destroy stability. This can be shown from the following estimate on the magnitudes of the RHS of the inequality for instability. We use the fact that i i w i i < \ / ^ ( i - * 2 r 1 / 2 - (A.210) V 7rn Appendix A. Bifurcations in an angular space 138 (See Abramowitz and Stegun (1970); Macrobert (1967).) then it follows that kn < -^L, kn{\ - kn) < -%=, (A.211) >n \/n where C is some positive constant of order 1. So, if n > nc, nc « A ll2 then the n ' th harmonic would not destabilize the homogeneous distribution. A . 3 Non- l inear analysis A .3 .1 Synerget ic analysis of M o d e l I in 2 D Under the assumptions given in section 2.6, we can perform the calculations shown below. First, form the expansion 'C{9,t)\ (C\ /£i(*)\ - (ii{t)\ = - + W*) + E W ) - (A-212) (Note that u;(#),i = 2 , . . . are cos(i9) for integer i, and (£;(£), r/,-(t)) are their amplitudes as functions of time, assumed to depend on the amplitude of the first mode. We consider the case when the initial angular distribution is even. The dynamics of the model preserves the evenness of the distribution.): m\ (umMt)]\ ,A , i = 2 , 3 , . . . (A.213) Mi)) WM*)^i(*)]/ A valid asymptotic expansion for {£,i{t), rn(t)) i = 2 , 3 , . . . in powers of (£1,771) is: 0 0 &(*)= £ cJitJa 6(*)Jl mi!)*. (A.214) Ji ,J2=0 A similar expression holds for r)i(t). I look for a solution of equations (2.43) in the form 0 0 C{0) = C + £{0), £(0) = £ ^ ( 0 ) , i eWI«C, (A.215) t = 0 Appendix A. Bifurcations in an angular space 139 i.e. the deviation away from the homogenous state is expressed as a superposition of the harmonics. Restricting attention to the linear approximation leads to equation (2.47). Substituting in the Fourier series expansion for £(#) given in equation (A.215), using the fact that K * vn = Knvn for every n = 1,2, . . . and equating coefficients of each harmonic on both sides of the resulting equation, results in the following system of equations: (-n2t+ aC-Kn{\ - Kn))in = 0 n = 1 , 2 , 3 . . . (A.216) (a - u) The expressions multiplying £„ in equation (2.47) coincides with the ones for the eigenval-ues of the stability matrix (or linear growth rates), so we can rewrite the above equation in the form: An£n = 0, \ n = -n2t+ aC Kn(l-Kn), n = l , 2 , . . . (A.217) [a — U) Note, that close to the primary bifurcation, Ai > 0, \ n < 0,n = 2, . . . . These equations have only trivial solutions, which indicates that the linear approximation is not sufficiently informative. Keeping terms up to second order in £ leads to a nonlinear equation. In this approximation equation (2.43) has the form: eA£ + a-f-[(K *t-K* (K * £)] + [^(K * £)2 + ft{K * fl - ft{K * (K * 0 ) -a-f-K * ((K * 02) - fK * U(K * 0)1 + [a^(K * 03 ~ J*1* * ((K * (f) -fK * U(K * 02) + f^K * ^ - f*K * W * ^ ~ 7*tK * &K * ^ = °' (A.218) where f = a — C. Substitute the expansion (A.215) into (A.218) and keep terms up to third power in £i and up to second power in £,-, i = 0 , 2 , . . . . We use trigonometric identities to simplify expressions involving products of the harmonics, (i.e. identities for products of sines and cosines) and equate coefficients of each harmonic on both sides of the equation as before. We get a hierarchical system of equations (2.48). Appendix A. Bifurcations in an angular space 140 The expressions for constants appearing in equations (2.48) and (2.49) are as follows: Al = ^ ( 1 - ATO, M = ^ ( 4 ^ ( 1 - ^ 0 - ^ ( 1 - ^ ) ) , B* = ^ Wl(2a - C)KI ~(a- C)Kf - CKl - (a - c)ki\, G = 2RF[2(2« - C)KX - (a - C)(^i + I<1) - 2CK, - (a - c)kx + ki)\, By = 2 ( ^ F [ 2 ( 2 « - c)kxk2 - (a - c){kxki + k2kD-2ckfk2 - (a - Qk^ki + k2)], r ~ 4(a-C)2\l ^ a-C A a-C >' D =GA-F + ^ , (A.219) where C = - ^ j . (A.220) Double h u m p e d kernel In the case of a double humped kernel the expansion is over the even harmonics only and has the form: oo C{9) = C + i\v2{0) + £ e2nv2n(0), (A.221) e2 = ±k'^\e-e>c\, | ^ | = ^ | c - < | B / 2 , n = 2 , 3 , . . . (A.222) where k',k'2n,e'c have the same meanings as before. The character of the bifurcations is the same. A.3 .2 Bifurcation Analysis of M ode l I in 3 D Consider an expansion of C in terms of the eigenfunctions (Legendre polynomials). oo C(ct>,0,t) = C + z0 + z^Pficoscj)) + z2{t)Pt{coscf>)cos{6) + J2yn(t)yn(<f>,e), (A.223) n=2 Appendix A. Bifurcations in an angular space 141 where the definition of the norm is: |Y„| = max|y n m | , (A.224) and where z\,z2 are amplitudes of the two unstable harmonics, and yn{t), ra = 2 , . . . are amplitudes of the stable harmonics. D o u b l e h u m p e d kernel In the case of the double humped kernel, the leading mode possesses both orientational and pattern degeneracy. Y2 = XlP° + x2P\ + x3P*. (A.225) The decomposition has the form: C(<f>, 9,t) =C + x0(t) + xx{t)P%{<£) + x2(t)P}(<t>) cos e + x3(t)P22(ct)) cos 20 + J2?n,n=2p2n(t)Y2n(t,e), PnYn = {p™Y™} m = 0 , . . . n , \p\=max\p™\. (A.227) As a result of intermode interaction, the pattern degeneracy is removed: all non-axisymmetric harmonics die out and if one choses again <f> = 0 as the direction of the axis of symmetry of the pattern evolved, then in the stable stationary state, x2 = x3 = 0, |p4J, \xo\ ~ x\, |p„| < x | , n > 4, and the leading mode amplitude obeys the equation Xxi + pxx\ - p2x\ = 0, (A.228) where A 2/^1 Pi = 4 9 ( a a _ g ) 3 ^ ( 1 - A '2) > ° ' ft ~ W > 0. (A.229) A . 4 Resu l t s for the Act in-Binding Kernel In this appendix I will briefly describe the results on orthogonal alignment. It is natural to assume, based on our experience with the single and double humped kernels, that Appendix A. Bifurcations in an angular space 142 an even, nonnegative double humped kernel having humps at —7r/2 and TT/2 leads to orthogonal alignment. (Alignment occurs along angles where the kernel has maxima). An example of this type of kernel is shown in Figure lc . Because of the periodicity of the kernel (period is 7r), only the coefficients Kn with n even are nonzero in the Fourier expansion for the kernel in 2D case, and expansion over spherical harmonics in the 3D case. It is important that both in 2D and in 3D, K2 < 0,K2(l — K2) < 0 and 0 < K4 < 1. This means that the second harmonic is stable and in general it will be the fourth harmonic that breaks stability and becomes the leading mode. (See also Civelekoglu and EK, 1994). This result holds for all three models considered in Chapter 2. In the 2D case, the bifurcational analysis close to criticality is analogous to the one carried out in Section 2.6. The bifurcation is supercritical, and C{0) * C + fc^/|e-ec|cos(40). See Figure 6c. In the 3D case, where again, the fourth harmonic Y4 is unstable, the leading mode is highly degenerate. Mode competition leads to removal of the pattern degeneracy, however, and the axisymmetric solution loses its preferential status. (It dies out as a result of competition.) The pattern evolved has the form of six mutually perpendicular smooth bumps in the angular distribution on the surface of the sphere. (See Figure 7c.) The bifurcation as in the 3D, double-humped case is transcritical character. It is natural to assume, based on the results of chapter 3, that far from criticality, where e —> 0 (diffusion is weak), four smooth bumps in the angular distribution generate four equal and mutually orthogonal narrow peaks, located at 6 = {0,±7r/2, n}. The shape of these peaks coincides with the shapes described in the case of the single humped kernel. Appendix A. Bifurcations in an angular space 143 In the 3D case, in the limiting case of weak diffusion, e —> 0, six equal, mutually orthogonal narrow peaks evolve from these bumps. I would conjecture that in the most convenient case, when one of the peaks is located at the pole, their location is given by the set of coordinates {</> = 0},{4> = 7r},{</> = ir/2;9 — 0,7r, — 7r/2,7r/2}. (The proof of this result is not complete.) Append ix B Calculations for the Peak Ansatz Derivat ion of equat ions (3.69) Plugging the expression P = (M — C) into the di-mensionless equation for C (the first equation in (2.20)) and assuming e = 0, we obtain (B.230) | 2 = _C[K * C) - C{K * M) + C{K * C) + a(M - C) = C(a + (K * M)) + aM. Plugging P = {M — C) into both equations (2.20) and adding leads to dM _ dP | dC dt dt "^ dt = P(K * C) - C(K * P) V ; V ' (B.231) = M(K *C)~ C(K * C) - C(K * M) + C(K * C) = M(K*C)-C{K*M). Proof that F(9) is constant The aim now is to show that F{6) is a constant, i.e., that F{6\) = ... = F(9n) where 9i...9n are the angles at which M ^ 0. I argue by contradiction: Suppose that F(9) is not constant. Then there is some angle 9\ such that Fx = F(91) > F(9), (B.232) (i.e., at this angle F attains a maximum.) for all 9 / 9\. Then from equation (3.75), it follows that i^ = E „a^ml\K^ ~e>)M{e^ (B'233) But on the right-hand side, the ratio a + Fl->l, 3 + l. (B.234) a + F{93) 144 Appendix B. Calculations for the Peak Ansatz 145 But then since K{6-6j) > 0 , M{93) > 0 , the right hand side of the above sum is greater than the expression n £ K{9X - 93)M(9J) = (K * M){9X) = Fx. 3=1 (B.235) (B.236) Thus we have Fx > Fi which is a contradiction, and therefore it must follow that F(9) is constant at the steady state. B . l Stabi l i ty of two peaks The linearized system of equations for the perturbations cx, c2, m has the form: m = 2KxCm - &&-*. + ^ c 2 , cx = (a + C(KX - K0))m - (a + ^f)cu (B.237) ( 62 = -(a + €{KX - KQ))m - (a + ^ ) c 2 , where C = Cx = Ci- We rewrite the system of equations (B.249) in the vector form where and where q = Aq . q = /m\ Cl , A = / & / w l-/ 6 =2A'1C, rf = ^ , / =(a + (7(A'1-e = -(a + £ 2 / • -d d\ e 0 0 e) K°))> (B.238) (B.239) (B.240) Appendix B. Calculations for the Peak Ansatz 146 Then a simple calculation reveals that the characteristic equation associated with the above matrix is det{A - XI} = (6 - A)(e - A)2 + df(e - A) + df(e - A) = (e - X)[{b - X)(e - X) + 2df] (B.241) = 0. Thus one of the eigenvalues, Ax is simply Ax = e = —(a + K^-) < 0 and the other two satisfy the quadratic equation [(b-X)(e-X) + 2df]=0. (B.242) An elementary calculation reveals that Xmax = ^(k -l2 + y/(k - Pf - 4kl'l), (B.243) where rs *, , KM ., {Kl-KQ)M ,nnAA, k = KxaM, l = a + ——, V = K 0J—. (B.244) Since A'i < KQ,1' < 0, and —kl'l > 0, it follows that at least one of the eigenvalues is positive. Thus, small deviations from the stationary steady state grow exponentially and two peaks are hence unstable. B.2 Stabil i ty of one peak I now consider the stability of one peak. The linearized system of equations for the perturbations m, ci,C2 has the form: m =K1€m + KiMc2, 6, =(a + C(K1-K0))m-(a + K0M)c1, (B.245) C2 = —am — (a + K\M)c2. Appendix B. Calculations for the Peak Ansatz 147 If we rewrite these in the vector form q = B q . where q is the column vector (m, c\, c-i) then the matrix B is (b 0 d\ B = e g 0 \f 0 */ Here b=K1C, d=K1M, e = a + C(Ki - K0), g = -(a + K0M), f = - a , t = -(a + KiM). We now find that the characteristic equation is det{B - AI} = (b - X)(g - X)(t - A) - df(g - A) = (g-X)[(b-\)(t-\)-df) = 0. (B.246) (B.247) (B.248) Thus Ai = g < 0 is one eigenvalue, and the others are solutions of the quadratic equation A2 - (b + t)\ + (bt - df) = 0. (B.249) This leads to A = i((6 + 0 ± J(b + ty + 4(df-bt)) = i [ - 7 ± yff + t], (B.250) where T<\ Kn M 2 n K, M2 (B.251) (a + A 'QM) (a + A'oM) Appendix B. Calculations for the Peak Ansatz 148 As 7 > 0, and £ < 0, both eigenvalues are negative. So one peak is stable to a small perturbation consisting of another competing peak. B.3 Stabil i ty of three peaks The masses at the three peaks satisfy the system of equations: Mx = Kx(Mi(C2 + C3) - Gx(M2 + M3)), M2 = K1(M2(Cl + C3)-C2(M1 + M3)), M3 =K1(M3(C1 + C2)-C3(M1 + M2)), Cx = - C i (a + A W i + KtM2 + Kx Ms) + aMx, C2 = -C2(a + K0M2 + KiMx + Kx M3) + aM2, C3 = -C3(a + K0Ms + A'iMi + K\M2) + aM3. A stability calculation in the case of these equations now leads to the 6 x 6 matrix (B.252) / where 26 - 6 -b d -b -b -b - 6 - 2 / / 26 - 6 / - 2 / - 6 26 / / - 6 - 6 - e 0 d - 6 0 - e - 6 d 0 0 6 = A^C, / = A'iM/3, d — a ~ KQC, e = a + A^M/3, / / - 2 / 0 0 —e \ (B.253) (B.254) Appendix B. Calculations for the Peak Ansatz 149 and where K = K0 + A'i + K\,C — J^ZWMT^- Observe that this matrix has a block structure. (a+KM/3)' 2b -b -b -b 26 -b -b -b 2b d -b -b - 2 / / / - 2 / / / or, -b d -b -b -b d /A f —e 0 0 B \ f - 2 / b o - e 0 (B.255) \C : D J where A, B, C, D are 3 x 3 matrices and D is diagonal. This fact is useful since it us to use a fact from matrix theory, namely that I A \ B\ (B.256) permits det det(AD - BC). (B.257) \C \ DI The eigenvalues of the block matrix can then be calculated by letting A' = A- M, D' = D- XI, in the matrix (A'D' — BC), with the result (B.258) (A'D' - BC) (X V fi\ V X V (B.259) where x = -(e + \)(2b-\)+2f(d + b), fi = (e + X)b-f(d + b). (B.260) Appendix B. Calculations for the Peak Ansatz 150 This leads to det(A'D' - BC) = (fi-X) [2^2 - l*X - X2} • (B.261) The equation 2fx2-nx-X2 = 0, (B.262) has four roots, and the case (fj, = x) gives two roots. An expanded form of this last equation is A2 + (e - 36)A + 3(/6 + fd - eb) = 0. (B.263) The roots of this equation are given by the expression: A = i ( - ( e - 36) ± y/(e - 36)2 + 4r , (B.264) where r = K\MC{KQ — K\) > 0. Thus at least one of the eigenvalues is positive, and three equal, evenly spaced peaks are unstable. Three unevenly spaced peaks As a specific example, we took the following arbitrary values, for the parameters: K(0) = 6, K{6X - 02) = 2, K(0! - 03) = 3, K(92 - 03) = 4, Mx = 10, M2 = 9, M3 = 4,. Then KM = 90. We assume a = 10, then Ci = l,C2= 0.9, C3 = 0.4. Using the software Mathematica, we found that in this case the 6 x 6 stability matrix has the form / 3 - 2 - 3 - 3 0 20 30 \ - 1 . 8 3.6 - 3 . 6 18 - 3 6 36 - 1 . 2 - 1 . 6 6.6 12 16 - 6 6 4 - 2 - 3 -100 0 0 - 1 . 8 4.6 - 3 . 6 0 -100 0 v -1 .2 - 1 . 6 7.6 0 0 100/ Appendix B. Calculations for the Peak Ansatz 151 Further, the six eigenvalues of this matrix are: Ai = -92.74, A2 = -96.92, A3 = 2 . 0 5 , (B.265) A4 = 0.00, A5 = 0.81, A6 = -100.0. Two of these eigenvalues are positive, so that this situation of three uneven peaks is unstable. Notice that one of the eigenvalues is zero since M = Mi + M2 + M3 = const. Append ix C Calculat ions for the spatio-angular case I consider Models I-III. Results are given below and summarized in Table I. M o d e l I Consider separately three cases, I: n = 0, II: n ^ 0, Kn > 1/2, III: n ^ 0, Kn < 1/2. Denote as g„(q) the right hand side of the inequality (4.144). Case I: The function -q2 -q2 9o{q) = e 2 (1 - e 2 ) has the asymptotes: 2 8 V + 0(<76), 9 « 1 <*>(<?)«< 2_2 8 (C.266) I e " ^ , <? > > 1 and the maximum value, go(q) = 1/4 at q = qo = (21n(2))1//2. If e2 > e2 = 1/2 then the LGR \o(q) < 0. If c2 < e2 then A0(g) < 0 whenever g > q^ and A0(g) > 0 whenever 0 < q < q^ where z(c) « ( i ^ ) 1 / 2 . (C.267) The left and right sides of the inequality (4.144) are shown on the LGR diagram (See Figure 17a). The LGR is maximal A-max ~ 7y ? (0.268) at q = qmax « q^/2 (See Figure 18a). Finally, we sketch the stability diagram in the governing parameter space (Figure 19a). Here the shaded region is the region of stability, (the other region corresponds to instability to the mode n = 0.) 152 Appendix C. Calculations for the spatio-angular case 153 C a s e I I : In this case, the LGR diagram is as shown in Figure 17b. The function gn(q) has the asymptotes 9 „ w 4 * ( 0 ) + A ; ( 2 A r 1 ) ^ ^ ¥ l l s l + o ( < ' 6 ) ' 9<<1 <c-^) [ -txn2 + K„e~g\ q » 1 gn(0) = - e m 2 + Kn(l - Kn) (C.270) and max<7„(g) = gn(0) + 1/4 at g0 = (21n(2^„))1 /2 . (C.271) If ei < e[c' where e^ is defined as - e[c)n2 + Kn{l - Kn) = 0, (C.272) then Xn(q) < 0 if q > q^ where q(c) ~ n (£ lZLl i l ) i /2 5 ( c . 273) at any e2. Otherwise, Xn(q) > 0. If ei > c[c' then Xn(q) < 0 at any q, whenever e2 > e 2 ( e i ) and Xn(q) > 0 at q m qo ^ 0 whenever e2 < e 2 ( e i ) 5 where ^ i ) * 1 " ^ 1 " 2 ' 4 (^ i ) < ^ n ( 2 A » - 1) < \ e.K—. (C.274) The LGR is maximal at q = 0 at ei < e[ (see Figure 18a). It is maximal if q « g0 7^ 0 if 1 4n2 ej < cj < — , e2 < e2(£i)- (C.275) (See Figure 18b). The stability diagram then has the form shown in Figure 19b. C a s e I I I : In this case, the LGR diagram looks like Figure 17c. Xn(q) > 0 at q < q^ whenever c\ < e\ and any e2. Xn{q) is the same as the one shown in Figure 18c. The stability diagram is shown in Figure 19c. Appendix C. Calculations for the spatio-angular case 154 To determine the primary bifurcation in the case that e\ and e2 decrease monotonously, we must find the intersection of a countable number of subsets of R2 corresponding to the shaded stability regions discussed above (for all n = 0,1,2...). Then the border of this intersection is the line of the primary bifurcation. The borders of the stability areas which do not belong to this line together form lines of secondary bifurcations. The cumulative set of bifurcation lines is shown in Figure 14. We call this picture the C o m p l e t e Bifurcation Diagram. M o d e l II It is easy to see that the right hand side of the inequality (4.146) has the following asymptotes: ((-ein2 + Cn2Wn) + C(Vn-iWn)q2 + 0(q% q « l 9n » < _ A _?2 (C.276) I -e^2 + CVnq2e-r-, q » \ Furthermore, the function gn(q) has a maximum at q = 0 if Vn < n2Wn/2 and at q ^ 0 if K > n2Wn/2. We consider separately three cases: (I) n = 0; (II) n ^ 0, Vn > n2Wn/2; (III) n^0,Vn<n2Wn/2; C a s e I: When n — 0 the LGR diagram is as shown in Figure 17a. Independently of ei there is stability at all q's if e2 > CVn. There is instability at small g's if e2 < CVn. The linear growth rate and the stability diagram are shown in Figures 18a and 19a respectively. C a s e I I : When n / 0, K > n2Wn/2, the LGR diagram is given in Figure 17b. If e2 > e2 a n d £i > ei \ where 4C> = C(Vn - ^ ) , 4C) = CWn, (C.277) there is stability at all g's. If e2 > e2 , t\ < e[ there is instability at small g's. At e2 < e2 , ei < ei < e i there is stability at all g's if e2 > e2(ei) and instability at Appendix C. Calculations for the spatio-angular case 155 q sa q(c) / 0 if e2 < e2(£i)- Here e^ and e'2(ei) depend parametrically on Vn, Wn, and C and can be found from certain transcendental equations, q^ is defined by the expression: q(c) ~ (2 _ ^ - f l \ (C.278) The linear growth rate and the stability diagram are given in Figures 17a, 18b and 19b. Case III: The LGR diagram in the case n ^ 0, Vn < n2Wn/2 is given in Figure 17c. Independently of £2, there is stability at all g's, if ei > e[c' and there is instability at small g's if ei < e[c' where - e[c) + CWn = 0. (C.279) (If Wn < 0 there is no instability at all.) The linear growth rate and the stability diagram are shown in Figures 18c, 19c respectively. Because e^ < CVn, the complete bifurcation diagram has the same form as the one in Figure 14. M o d e l III The right hand side of the inequality (4.148) has the asymptotes: / ( - £ l n 2 + r)C\l - Gn)) + \C>Gnt + 0{q% q « l 9n(q) ~ < _ _ A a2 (C.280) \{-t1n2^r]C2)-riC2Gne-ar. q » \ gn(q) is monotonously decreasing if Gn < 0 and monotonously increasing if Gn > 0. Here we have to distinguish three cases: (I) n = 0; (II) n ^ 0, Gn > 0; (III) n / 0, Gn < 0. Case I: The LGR diagram for the case n = 0 is shown in Figure 17a. Independently of ei there is stability at all g's if t-i > r/C2/2 and there is instability at small g's if e2 < r]C2/2. The linear growth rate and stability diagram are given in Figures 18a, 19a respectively. Case II: If n ^ 0, Gn > 0, the LGR diagram is similar to that shown in Figure 17b. The only difference is that the humped curves are monotonically increasing in this case, Appendix C. Calculations for the spatio-angular case 156 a feature which does not alter the results. If £2 > 4 a n d £1 > 4 , where 4=' = ^ , £ l(c) = „ ( - < i ^ i , (C.281) there is stability at all g's. At £2 > 4 and t\ < e[ , there is instability at small g's. At £2 < 4 1 ei < £i < f]C2/n2 there is stability at all g's. If £2 > £2(ei) there is instability at some o ~ q(c) =jL 0 if £2 < £2(ei)- Here £2(ei) depends parametrically on rj, C, Gn, and can be found from certain transcendental equations. The distinctive feature of this model is that q(c> can vary from zero to infinity: nC 2 q(c> -> 00 if £2 ->• 0, £1 -+ - L T - . (C.282) n2 The linear growth rate and stability diagram are shown in Figures 18a,18b,19b. C a s e I I I : If n 7^ 0, Gn < 0, the LGR diagram is shown in Figure 17c. Independently of £2, there is stability at all q's if £1 > 4 and instability at small q:s if e\ < 4 where 4 is given in (C.293). The linear growth rate and stability diagram are shown in Figures 18c and 19c, respectively. Because 4C) = ^ G n < *£-, (C.283) the complete bifurcation diagram has the same form as the one in Figure 14. Append ix D Supplement to the linear stabil ity analysis of Chapter 5 D . l Propert ies of the bifurcation function (5.176) In this appendix we describe briefly the properties of the bifurcation function defined by the expression (5.176): %) = -rrw - -rr-2- (D-2 8 4) qz + bz q2 -f a2 These are based on elementary calculus. The bifurcation function is even, and thus considered at positive values of the argument. Regions 1 and 10 (B < A:a(B/Ay'4 < b): The bifurcation function is negative and increasing and has an absolute minimum at the origin (Figure 23b). Region 2 (B < A,a(B/A)1/2 <b< a{B/Afl4): The bifurcat ion function is negative and has an absolute minimum at q* = (B*a2 - Ah2)* (A? - B*)~* . (D.285) (See Figure 23f.) Regions 3 and 4 (B < A,b < a(BfA)1/2): The bifurcation function is positive and decreasing at small values of the argument, and negative at large values of the argument. It has a positive maximum at the origin and a negative minimum at qc given by (D.285) (Figure 23c). Regions 5 and 6 (B > A,b < a(B/A)1/4): The bifurcation function is positive and decreasing and has an absolute maximum at the origin (Figure 23a). 157 Appendix D. Supplement to the linear stability analysis of Chapter 5 158 Region 7 (B > A,a{B/A)llA < b < a(B/A)1/2): The bifurcation function is positive and has an absolute maximum at qc given by (D.285) (Figure 23e). Regions 8 and 9 (B > A,b > a(B/A)1'2): The bifurcation function is negative and increasing at small values of the argument, and positive at large values of the argument. It has a negative minimum at the origin and a positive maximum at qc given by (D.285) (Figure 23d). D.2 Propert ies of the LHS of the instabil i ty criterion (5.188) The primary bifurcation in the "finite-range" case corresponds to a positive maximum of the LHS W2(q2)~X(qi), (D.286) of the inequality (5.188). Our goal is to find the critical wave vector qc = (qic,q2c) corresponding to this maximum. Regions 1, 10 (B < A,a(B/Ay/4 < b): The bifurcation function is negative and has an absolute minimum at the origin {q\c — 0). Thus, the maximum of expression (D.286) corresponds to the negative minimum of the function W2(q2), or to q2c — 4.5. Region 9 (B > A, a(B/A) < b): Expression (D.286) has two local maxima correspond-ing to qx = 0,<?2 - 4.5 and qx = qc*,q2 = 0. The inequality J^2(4.5)A(0) > W2(0)X(qc*) is easily established giving qic = 0,q2c — 4.5. Region 8 (B > A^a(BfA)1^2 < b < a(B/A)): The bifurcation function is the same as in Region 9. However, if b > WA,a(B), where b = WA,a(B) is the equation of the curve (dashed line in Figure 22) defined by the equation: ^2(4.5)A(0) = W2(0)\(qc*), (D.287) Appendix D. Supplement to the linear stability analysis of Chapter 5 159 then the situation is equivalent to the one in the previous paragraph. Equation (D.287) is equivalent to the equation (5.189): recall that W2{q2c) = W^(4.5), \{qic) = X(qc*) and W2(0) = 1. If b < wA,a(B) then T^2(4.5)A(0) < W2(0)X(qe*), and qlc = q*,q2c = 0. Region 7 (B > A,a(B/A)1/4 < b < a(B/A)1/2): The bifurcation function is positive and has a positive maximum, thus the maximum of the expression (D.286) corresponds to qlc = qc*,q2c = 0. Regions 5 and 6 (B > A,b < a(JB/A)1/4): The bifurcat ion function is positive and has an absolute maximum at the origin, so the maximum of the expression (D.286) corresponds to q\c = 0, q2c = 0. Region 4 (B < A,b < a(B/A)): The expression (D.286) has two local maxima cor-responding to q1 = 0, q2 = 0 and qi = q*-,q2 — 4.5. The inequality 1^2(0)^(0) > W2(4.5)\(qc*) is easily established giving qic = 0, q2c = 0. Region 3 (B < A,a(B/A) <b< a{B/A)ll2): The bifurcation funct ion is the same as in Region 4, however, if b < WA,a(B), then the situation is equivalent to the one in the previous paragraph. If b > WA,a(B) then l^2(0)A(0) < W2(4.5)X(qc*), and qic = qc*,q2c = 4.5. Region 2 (B < A,a(B/A)1/2 < b < a(B/A)1/4): The bifurcation function is negative and has an absolute minimum at q\c — q*. Then expression (D.286) has a maximum at <?ic = qc*, q\2c - 4 . 5 . Table I Angular mode Model Case I n = 0 Case II I: Kn > 1/2 ILK, > n2Wn/2 III:G„ > 0 Critical wavenumber f i £2 Scenario Figure Growth rate diagram (Fig) Linear growth rate (Fig) Stability diagram (Fig) q(c). any (2 < B 15b 17a 18a 19a ~ 0 c(<0 e2 7<c) = 0 e2 > e<e) A 15a 17b 18a 19b «<c> £ 0 M) < ei < ei e2 < e<e) C 15c 17b 18b 19b Case III l:Kn < 1/2 Il:V„ < RHS III:G„ < 0 ] ( c ) = 0 *i < e « any A 15a 17c 18c 19c 160 Interaction kernel (Fig) 23c 23c 23c 23c 23b,f 23d 23d 23d 23d 23a Bifurcation function (Fig) 23b 23f 23c 23c 23a 23a 23e 23d 23d 23b ID or 2D "con." case: critical wave vector; form of pat tern no pat tern no pat tern qic ~ 0 , g 2 c = long-ranged qic ^0,q2c = long-ranged qic ^ 0,g2c = long-ranged qic ^ 0,q2c = long-ranged qic ^0,q2c = rolls qic ^ 0,g2 c = rolls 91c 7^0, q2c = rolls no pat tern 0 • 0 0 0 0 0 0 2D"finite-ranged' case: critical wave vector; form of pat tern qu = 0,q2c # 0 pat tern - ? qic ^ 0 , g 2 c ^ o rippling qic ^ 0 ,g 2 c ^ 0 rippling qic ~ 0,g2 c = 0 long-ranged qic ~ 0,92c = 0 long-ranged qic ~ 0,92c = 0 long-ranged qu ^ 0,92c = 0 rolls qu ^ 0,92C = 0 rolls 9ic = 0,92c ^ 0 pat tern - ? 9ic = 0,92c ^ 0 pa t tern - ? 161 (a) -Jt - n / 2 n/2 7t (b) -re - r t /2 0 TC/2 (C) -7C - r t /2 0 7C/2 Fig. 1 The kernels K used to represent the contact alignment phenomena are shown here as functions of the angle between cells. (In 2D the angle is 8, and in 3D it is 7.) (a) Single humped kernel, representing alignment in which both cells are oriented head-to-head, (b) Double humped kernel which permits cells to align also m a head-to-tail configuration, (c) An orthogonal interaction kernel that is not relevant for cellular interactions but that plays a role in the model for actin alignment. 142, Fig.2 (a) The expression K„{\ — A„ ) which appears in the condition for instability for Model I is shown plotted as a function of the mode n for the kernei A'/ Note that though the curves are shown as if they were continuous we are only interested in behaviour at integer values of n. The figure is generated using Mathematica to calculate the Legendre coefficients of the kernel. Note that only modes n — 1. 2. 3 can rnuse instability (b) The same expression, but for the kernel A"//, (c) for Km- We also show the left-hand side of the inequality, i.e.. An{rt -<- 1} as a function of n for A — 0.1 and A = 0.03. Note lhai m this case the1 modes n — 1 and » = 2 will satisfy the condition for instability. The critical value of .4. i e the one causing the onset of instability i- A„ ~ 0.0S in (a), ACT a 0.06 in (b). and ,4 , r % 0.04 in (c). 163 Fig.3 (a) The expression CWn which appears in the condition for instability for Model II is shown plotted as a function of the mode n for the kernel I in 2D. (b) The same expression as in Figure 3a. but shown together with superimposed lines representing two possible values of £, namely e = 1.0,0.7. It can be seen that the inequality criterion is satisfied in the latter case and instability to the mode ;j = 1 occurs. (c) As in 3(a) but for the kernel I in 3D. (d) Same as 3(b) but in 3D. The values <= = 0.3,0.2 are shown. Instability to n = 1 occurs for e = 0.2. U4-Fig.4 (a) The expression rjC1^ — G„) which appears in the instability condition for Model III is shown plotted as a function of the mode n for the kernel I in 2D. (b) The same expression as in 4a, with a superimposed graph of the expression en2 for two possible values of e, namely c = 0.2, 0.15. The smaller e value causes instability to n = 1. (c) As in 4a, but for the kernel II in 3D. (d) The expression in 4c. We also show the expression en(n + 1) for two possible values of e, namely c = 0.2,0.4. 165 (a) * i (b) • 8 Fig.5 (a) Supercritical bifurcations found in the analysis of the 2D and 3D cases, (b) Transcritical bifurcation in the 3D double-humped kernel case. 166 (a) (b) % 2n (c) n 2n Fig.6 Types of patterns seen in the cell distribution in the 2D case after bifurcation, (a) A distribniion similar to cos 9 which occurs in the case of a single-humped kernel, (b) A distribution similar to cos 19 occurs when the kernel is double-humped, (c) A distribution similar to cos 49 occurs in the case of orthogonal alignment kernels. 167 (a) (b) (c) 1'ig.T Three surface spherical harmonics that are significant in growing patterns in the 3D models, (a) P-'. (b) f2°. (c) Y'43: arises in the actin binding model where orthogonal interactions occur. The shaded regions on the sphere represent locations of increased density. 168 (a) Fig.8 We simulate the evolution of peaks for Model I. Development of an initial distribution (dashed lines ) consisting of three small peaks under the dynamics predicted by the model (2.1). Shown is the density of bound cells only. (The density of free cells is quite similar, but of different magnitude.) The turning due to cell contact is here based on a strictly positive, single humped kernel K{6) - 1.5 + cos 9. The discretised equations were simulated for ,9 = 7 = 0.3, e = 0.001, At = 0.005 for (a) the first 3000 iterations, (b) the next 3000 iterations. A single peak forms, and stabilizes. 169 Fig.9 As in Fig 8 but with a single-humped kernel of the form (3.116) which is nonzero only on a compact set. Two peaks persist. All parameters and conditions were the same as in Figure 8. 170 180 -theto. 360 Fig 10 As in Figure 9, but with a kernel whose range is 55°. (See text for the exact form.) Note that now all three initial peaks can grow and coexist, since their interaction (mainly through angular diffusion) is very weak. 171 180 the-ta 360 Fig, 11 Same as Fig 10, but with an increased rate of diffusion, and simulated for longer time, t - 1.0, 15,000 iterations. Increased diffusion enhances the interactions of the peaks, and leads to the above situation. 172 (a) (b) Fig.12 We used a doubie-humped kernel K(S) = 1.5 +cos 20, and difi'usion rate t = 0.001. (a) after 3000 iterations, (b) the next 3000 iterations. Two unequal peak.s evolved. Because the difi'usion is so small, the peaks do not communicate with one another. 173 180 -theta 360 Fig.13 Similar to Figure 12, but with large diffusion, t = 1. Two peaks at relative angle rr evolve. Their heights are gradually equalized by the rotational diffusion. (Here shown up to 15,000 iterations). 174 k I§l Fig. 14 A complete bifurcation diagram. A and B are the lines of primary bifurcations. Line A corre-sponds to the growing mode n ^ 0,q = 0. Line B corresponds to n = 0,</ ~ 0. The shaded region is the region of stability of the homogeneous state. The unshaded region corresponds to some non-homogeneous pattern. The dashed lines are secondary bifurcations. The vertical dashed lines correspond to n / 0,<] = 0. The curved portion of the dashed lines corresponds to n / 0 , i ; ^ 0. 175 7/ /// V UJUJUUU \\\\ww\\\ tttttttttttt Fig.15 Depending on the growth protocol, any one of three possible bifurcation scenarios can take place Star t ing from an initially disordered state (left), bifurcation may lead to (A) formation of angular order in a distribution that is spatially homogeneous, (B) aggregation and formation of spatial inhomogeneities without angular order taking place, (C) formation of patches of aligned objects. 176 ^ Fig.16 The same as Figure 15, but with an additional line for primary bifurcation (C). This line corre-sponds to the growing mode n / 0,</ ^ 0. 177 (a) (b) * q El l a I ? e RHS (c) LHS RHS Fig, 17 The linear growth rate diagrams. The dashed lines show the left-hand sides of the instabili ty criteria of the three models. Solid lines are the right hand sides of these criteria. Instability occurs whenever the dashed line is below the solid line, (a) Case I: the unstable mode is n = 0, q = c/(c< « 0. (b) Case II: depending on the values of t\,t2, two possibilities occur. If £2 is large, and ei is small, we have a s i tuat ion similar to (a) above. If <r2 is small and ex is large, the mode n jt 0, q = ?(c> ^ 0 is the leading one. (c) Case III: the mode n ^ 0, q - 0 is unstable. 178 (a) Mq) Kti) (b) • q Uq) (c) > q Fig. 18 The linear growth rate of an unstable mode as a function of q. (a) and (<:) correspond to 17(a).(c) respectively, (b) corresponds to 17(a) if c2 is large, and to 17(b) if f2 is small 179 (a) • £ i (b) (c) _ _ >(cj • > £ (C) : - ^ - L - s - ^ - ; - . -A- *.,.,Tr * • • e (c) Fig. 19 The stability diagram in the governing parameter space. The shaded region is the stability region and the unshaded region corresponds to a non- homogeneous pattern, (a) Case I: Independently of ci, when (2 decreases below e[c\ the orientationally disordered patches start to grow (scenario B). (b) Case II: the vertical border of the shaded region corresponds to the growing mode n ^ 0,q = 0 (scenario A). The curved border corresponds to the growing mode n j ^ O , } / 0 (scenario C). (c) Case III: Independently of f2, as f i decreases below q . the mode n ^ 0,q = 0 breaks the stability (scenario A). 180 Fig.20. The velocity v of the given cell at the point (x,y) is caused by the interaction with the cell at the point (r'.y1). If the latter is in the domain 1, the interaction is repulsive (or attraclive, depending on the model parameters). If it is in the domain 2, the interaction is attractive (repulsive). In the unshaded region the interaction is exponentially small. 181 — I-Fig.21 A sketch of the travelling wave pattern of the bacterial swarm. 182 Fig.22 The parameter space (range of the attraction on y axis; strength of attraction on x axis) The space is divided into 10 segments. The value of the interaction parameters from the different segments correspond to the qualitatively different forms of the interaction kernel, bifurcation function and individual potential function. Explanations in the text. 183 K.X K.X x.q " < • • • < / K.X K.X x,q K, X t K.X •t.q Fig.23 The qualitative forms of the functions A"(.r) and X(q) is shown at positive values of the arguments for the interaction parameters in various regions in the parameter space. Explanations in the text. 184 d Fig.24 The periodic patterns of stripes and spots forming in two dimensional bacterial swarm. The pattern on Fig.24b resembles the rippling structure.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Modelling spatio-angular patterns in cell biology
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Modelling spatio-angular patterns in cell biology Mogilner, Alexander 1995
pdf
Page Metadata
Item Metadata
Title | Modelling spatio-angular patterns in cell biology |
Creator |
Mogilner, Alexander |
Date Issued | 1995 |
Description | In this thesis I investigate the formation of spatio-angular order in populations of interacting individuals. Mathematical models are used to describe distributions of orientation, and how these evolve under different assumptions about the turning behaviour of individuals. Three types of models are introduced: in one, turning is described by a gradual direction shift (rather than abrupt transition). In two other variants, interactions between individuals cause abrupt turning and mutual displacements. Linear stability analysis and weakly non-linear (synergetics) analysis of interacting modes are used to understand the nature and stability properties of the steady state solutions. I investigate how nonhomogeneous pattern evolves close to the bifurcation point and find that individuals tend to cluster together in one direction of alignment. A special limiting case of nearly complete alignment occurs when the rotational diffusion of individual objects becomes very slow. In this case, motion of the objects is essentially deterministic, and individuals tend to gather in clusters at various orientations. To understand analytically the behaviour of the deterministic models, I represent the angular distribution as a number of 8-\ike peaks. For weak but nonzero angular diffusion, the peaks are smoothed out. The analysis of this case leads to the study of a singular perturbation problem. Next, I generalize the models to examine the dynamic behaviour of ensembles of individuals whose interactions depend on angular orientations as well as spatial positions. I show how processes such as mutual alignment, pattern formation, and aggregation can be described by sets of partial differential equations containing convolution terms. Such models appear to contain a rich diversity of possible behaviour and dynamics, depending on the details of the convolution kernels involved. The analysis of the equations, and predictions in several test cases are presented. Finally, I introduce a partial integro-differential model of aggregation phenomena. The model can be derived from a system of chemotactic partial differential equations (PDE). A linear stability analysis reveals that under certain conditions, periodic patterns may emerge. A weakly non-linear analysis predicts stripe patterns. In the limit of slow diffusion I use a Lagrangian approach to derive and analyse the gradient system and obtain periodic peak-like patterns. The possibility of existence of other types of patterns is discussed. I consider the application of the model to a phenomenon known as rippling in a swarm of social bacteria observed in Myxobacteria. |
Extent | 7173179 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-04-24 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0079983 |
URI | http://hdl.handle.net/2429/7546 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1995-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1995-060268.pdf [ 6.84MB ]
- Metadata
- JSON: 831-1.0079983.json
- JSON-LD: 831-1.0079983-ld.json
- RDF/XML (Pretty): 831-1.0079983-rdf.xml
- RDF/JSON: 831-1.0079983-rdf.json
- Turtle: 831-1.0079983-turtle.txt
- N-Triples: 831-1.0079983-rdf-ntriples.txt
- Original Record: 831-1.0079983-source.json
- Full Text
- 831-1.0079983-fulltext.txt
- Citation
- 831-1.0079983.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079983/manifest