text:
# The Evolution of Accreting Black Holes in Outburst

## 1 Objective

The primary objective of this presentation is to describe our current understanding of the observational properties of black hole X-ray spectra. In addition to describing the phenomenology of black hole spectral states, I discuss physical implications of the observations, including how they may constrain accretion geometries and emission mechanisms. The flexible scheduling and broad bandpass of RXTE have been critical to these studies and especially to our group (Tomsick, Kalemci, Corbel, and Kaaret) as we have focused on frequent RXTE and radio observations close to state transitions. Much of the recent debate on black holes has been related to the physics of the Hard State, and, here, I provide a summary of how the predictions of the current theoretical models compare to observations. This presentation includes some discussion of black hole X-ray timing and radio emission, but these topics are dealt with in more detail by others in these proceedings (see papers by Remillard and Corbel). Finally, I discuss the instrumental capabilities that will be important in future studies of black hole spectra.

## 2 Black Hole Systems Observed by RXTE

During the RXTE lifetime (1996-present), 28 accreting black hole or black hole candidate systems have been observed in outburst. Compact object mass measurements were obtained, via optical or infrared observations, for 11 of these systems, and the masses are higher than 3$`\mathrm{M}_{\mathrm{}}`$, the theoretical mass upper limit for neutron stars, in all of these. Eight of the 28 sources have been in outburst for the entire RXTE mission, and Figure 1 shows an example light curve for one of the “persistent” sources (Cyg X-1). Four of the sources have had multiple outbursts between 1996 and now, and the light curve for the most active “recurrent transient” (4U 1630–47) is shown in Figure 1. The other 16 sources have had only one outburst observed by RXTE, but it should be noted that several of these sources had outbursts before 1996 (Chen et al., 1997) and are know to be recurrent.

## 3 Spectral States

Although the energy spectra of accreting black holes are complex, they are dominated by two emission components: a soft thermal component and a hard power-law or cutoff power-law component extending to hundreds of keV. The soft component is very likely blackbody emission from an accretion disk, and the disk may have properties similar to (or the same as) a Shakura and Sunyaev (1973) disk. Spectral states are defined, in large part, based on the relative strength of the two components. Recently, the large number of RXTE observations of accreting black holes in various spectral states were used to quantitatively define the spectral states (McClintock and Remillard, 2003). In addition to giving the most precise spectral state definitions to date, McClintock and Remillard (2003) considered the fact that there are many examples (see below) for which the traditional luminosity-based state names are not appropriate and re-named the spectral states. A listing of the three outburst states and their general properties follow (see McClintock and Remillard, 2003, for the complete quantitative definitions of the states). Figure 2 shows examples of spectra for sources in these three states.
$``$ The Steep Power-Law (SPL) state was formerly known as the Very High state. Both disk and power-law components are usually present in this state, and a source is said to be in the SPL state if either the power-law accounts for more than 50% of the 2-20 keV flux or if quasi-periodic oscillations (QPOs) are present and the power-law accounts for more than 20% of the 2-20 keV flux. In this state, the power-law component has a photon index of $`\mathrm{\Gamma }>2.4`$.
$``$ The Thermal-Dominant (TD) state was formerly known as the High-Soft state. The disk component dominates the X-ray spectrum, accounting for more than 75% of the 2-20 keV flux. The level of timing noise in this state is usually very low.
$``$ The Hard state was formerly known as Low-Hard state111Even McClintock & Remillard (2003) continued to use the old terminology for this state, but several examples indicate that the “Low-” should be dropped.. The power-law, with $`1.5<\mathrm{\Gamma }<2.1`$, dominates in this state, accounting for more than 80% of the 2-20 keV flux. There is a high level of timing noise, and radio emission, likely from a compact jet, is typical of this state (Fender, 2001).

An “Intermediate State” has also been defined as one of the canonical black hole spectral states (Mendez et al., 1998); however, the only property that clearly separates this state from the SPL state is that the Intermediate State has traditionally been used to describe black holes at lower luminosity. Thus, the argument that the states are not strictly dependent on luminosity makes this state, as canonically defined, unnecessary. However, there are times, especially during state transitions, where the source properties do not fit into the SPL, TD, or Hard states, and I agree with the suggestion of McClintock and Remillard (2003) that the states that show various combinations of SPL, TD, and Hard state properties be referred to as intermediate states.

There is a final state that is simply defined by (low) luminosity called the “quiescent” state. For black hole systems, a source is typically said to be in quiescence if its X-ray luminosity is in the $`10^{3033}`$ erg s-1 range (Garcia et al., 2001; Hameury et al., 2003; Tomsick et al., 2003a). However, due to the relatively poor statistical quality of the low luminosity observations, it is difficult to determine if there are other source properties that distinguish quiescent properties from Hard state properties, and, below, I discuss this issue further.

## 4 Examples of Luminosity Independence of States

Here, I provide some examples of sources that exhibit X-ray properties that fit in well with the standard spectral state definitions outlined above, but do not show the luminosity dependence that is considered to be typical. This is not meant to be an exhaustive list, but I simply discuss some of the clearest examples. It has long been known that some sources exhibit “Hard state outbursts” where they reach high luminosities, but remain in the Hard state. The most dramatic example of this behavior is seen from V404 Cyg (= GS 2023+338), which remained in the Hard state while reaching a luminosity of $``$$`2\times 10^{39}`$ erg s-1, assuming a distance of 3.5 kpc (Tanaka and Shibazaki, 1996). A more recent example is the behavior of XTE J1650–500 as this source began its outburst in the Hard state and remained in the hard state for two weeks (see Figure 3), during which time the 3-20 keV luminosity peaked. The source eventually made a transition to the SPL and TD states (Homan et al., 2003), and, later in the outburst, the source showed a transition back to the Hard state (Kalemci et al., 2003; Homan et al., 2003). Chandra observations occurred near the end of the outburst, showing that the source exhibited Hard state properties at a luminosity of $`9\times 10^{34}`$ erg s-1 (Tomsick et al., 2003b). Overall, XTE J1650–500 exhibited Hard state properties over a range of three orders of magnitude in X-ray luminosity. Another excellent example of luminosity independence comes from a study of a full outburst from XTE J1550–564 (Homan et al., 2001). The hardness-intensity diagram for XTE J1550–564 indicates that this source can enter the Hard state, the SPL state, or various intermediate states over a large intensity range. Other examples include GRS 1758–258 and 1E 1740.7–2942, which show transitions from Hard to TD states as the source luminosity drops (Smith et al., 2001, 2002). Also, Cygnus X-1 Zhang et al. (1997) and SAX J1711.6–3808 (Wijnands and Miller, 2002) illustrate the luminosity independence of spectral states.

In addition to being important for determining the most meaningful classification of spectral states, this luminosity independence implies that the states are not simply set by the mass accretion rate. As discussed by Homan et al. (2001) and others, this indicates that another physical parameter that is not directly tied to the mass accretion rate must be important in setting the emission properties. Although there have been suggestions for what this second parameter might be, such as the inner radius of the optically thick accretion disk ($`R_{\mathrm{in}}`$), the size of the corona, or the jet power, it is currently unclear which of these parameters cause significant changes in the X-ray emission properties.

## 5 X-Ray Properties and Physical Implications

In this section, I describe the specific observational properties of different phases of the evolution of accreting black holes and discuss the physical implications of the observations. The different evolutionary phases include: states that typically exhibit a significant or dominant soft component (SPL and TD); transitions to the Hard state; and the Hard state.

### 5.1 Steep Power-Law and Thermal-Dominant States

The SPL and TD state spectra both typically include a significant contribution from an optically thick disk, but, with the stronger power-law component in the SPL, a third “reflection” component arises due to the hard flux illuminating the accretion disk (Lightman and White, 1988). The two main features of the reflection component are an iron K$`\alpha `$ emission line near 6.4 keV and a “reflection bump” between 10 and 100 keV, and, for many broadband black hole spectra, these features are present (e.g., Gierliński et al., 1999). Recent observations by Chandra and XMM-Newton (e.g., Miller et al., 2002) indicate that many black hole systems have an intrinsically broad iron line, and the broadening is interpreted as being caused by a redshift due to the black hole’s gravitational field. The inner radii inferred from spectral modeling of the broad lines provides sone of the strongest evidence that the inner radius of the disk reaches the Innermost Stable Circular Orbit (ISCO). Other evidence that the disk extends close to the black hole comes from the detection of high frequency quasi-periodic oscillations (QPOs), which are only seen in the SPL state (see the paper by Remillard in these proceedings), and also from modeling the soft spectral component using a Shakura and Sunyaev (1973) model (although it has been shown that, in some cases, the inner radii inferred from such fits give unphysically small values).

While there is a relatively good understanding of the origin of the soft component and of the geometry of the optically thick portion of the accretion flow, the origin of the power-law component is currently unclear. However, the fact that the power-law in the SPL state can extend up to energies approaching at least 1 MeV (e.g., Tomsick et al., 1999) indicates that a non-thermal electron energy distribution is required. Understanding the production of this component in detail is an important, unsolved problem (see the paper by Coppi in these proceedings for more on this issue).

### 5.2 Transitions to the Hard State

A main focus of our group has been to measure the properties of accreting black holes in X-rays with RXTE and in the radio band close to state transitions. Frequent (e.g., daily) monitoring is critical for this study because it is known that state transitions can occur on time scales of days (Tomsick and Kaaret, 2000). The most accessible and reliable state transition is the transition from the TD or SPL state to the Hard state at the ends of transient outbursts, and performing radio observations close to such a transition is especially interesting for understanding the compact jet turn-on.

Using our observations as well as observations from the RXTE archive, Kalemci (2002) analyzed data from observations of 16 outbursts from 11 different sources, and, in eight cases, good RXTE coverage was obtained close to the state transition. The data are shown in the paper by Kalemci et al. in these proceedings, and I summarize some of the basic results here. The sharpest changes were seen in the timing properties, with sources typically showing a change from a very low, undetectable, level of timing noise to an RMS level of tens of % in less than a day. It is clear that the timing noise can turn on very rapidly in these systems. While the timing changes are very sharp, gradual changes are seen in the spectral parameters, including the power-law photon index, $`\mathrm{\Gamma }`$, and the temperature at the inner edge of the accretion disk, $`kT_{\mathrm{in}}`$. Finally, a very large drop in the 3-25 keV flux of the soft component is apparent, and, typically, this component becomes undetectable after the transition.

The sharp timing transition that is seen provides a clearly defined transition time, and Figure 4 shows radio measurements that were made close to this time for four systems. For 4U 1543–47, XTE J1550–564, and XTE J1650–500, we obtained radio detections after the state transition. For 4U 1630–47, we obtained an upper limit very close to the time of the state transition. The most interesting case is 4U 1543–47 because, when our data are combined with the measurements reported by Park et al. (2003), very good radio coverage was obtained. For this source, the radio emission increased, but not until 8-12 days after the state transition (see also the paper by Buxton et al. in these proceedings for information on the behavior of 4U 1543-47 during this time in the infrared). Our analysis of the RXTE data for 4U 1543–47 indicates that the first radio detection occurred about the time that the source reached the true Hard state as indicated by $`\mathrm{\Gamma }`$ reaching its hardest value of about 1.6 (Kalemci et al., in preparation). This suggests that the processes that lead to an increase in timing noise do not automatically result in jet production, and re-enforces the connection between jets and the hardness of the X-ray spectrum (e.g., Corbel et al., 2000).

### 5.3 Hard State

A current topic of debate concerns the questions of what is the X-ray emission mechanism and what is the source geometry in the Hard state. The standard picture is that the cutoff power-law is produced by inverse Comptonization from a thermal distribution of coronal electrons where the coronal temperature is close to the cutoff energy. While the energy spectrum is well-described by thermal Comptonization, the discovery that a compact jet is present in the Hard state has led to the question of whether at least some of the X-ray emission could be produced in the jet, and jet-based synchrotron (Markoff et al., 2001) and Compton (Georganopoulos et al., 2002) models have been investigated.

One specific issue concerning the source geometry is whether the inner radius of the optically thick accretion disk increases in the Hard state. One might think that the drop in flux and temperature of the soft component during and after the transition to the Hard state described above demonstrates that $`R_{\mathrm{in}}`$ increases; however, other possibilities have been suggested such as the idea that the inner disk remains intact, but, in the Hard state, the viscosity mechanism in the inner disk turns off and the accretion energy goes into producing the jet or heating the corona rather than being radiated away in the inner disk (Livio et al., 2003). Still, even if the inner disk does not radiate efficiently, one might expect to see a reflection component if the inner disk remains at the ISCO. Studies of the strength of the reflection component have been carried out, showing a large decrease in the strength of the reflection component as $`\mathrm{\Gamma }`$ becomes harder when the source passes through intermediate states on its way to the hard state (Zdziarski et al., 1999, 2003). While the trend in the strength of the reflection component is consistent with an increase in $`R_{\mathrm{in}}`$, other possible explanations for the trend include ionization of the inner part of the disk (Nayakshin, 2000), but also see (Barrio et al., 2003), or beaming of emission away from the disk (Beloborodov, 1999).

While the evolution of the soft component and the reflection component do not provide a definitive answer to the question of how $`R_{\mathrm{in}}`$ evolves, it is notable that changes in both of these components are well-explained by an increasing $`R_{\mathrm{in}}`$. Furthermore, an increasing $`R_{\mathrm{in}}`$ can provide a rather natural explanation for observed changes in timing properties. As shown in Figure 5, the black hole sources appear to show a uniform drop in QPO frequency as $`\mathrm{\Gamma }`$ hardens (Kalemci, 2002), and similar trends have been reported for either QPO frequencies or break frequencies in the power spectrum in other work (di Matteo and Psaltis, 1999; Nowak et al., 2002; Revnivtsev et al., 2001; Vignarca et al., 2003). An increasing $`R_{\mathrm{in}}`$ may help to explain these trends as it provides an increasing dynamical time scale in the system.

Recently, we used RXTE and Chandra observations of XTE J1650–500 to show that the trend of characteristic frequencies dropping with luminosity extends to very low luminosities near $`1\times 10^{34}`$ erg s-1 as shown in Figure 6 (Tomsick et al., 2003b). We detect a change in the break frequency by a factor of $``$1200 between the transition to the Hard state and the lowest luminosity we sampled, and if this is set by a dynamical time scale that is determined by the inner radius of the disk, the implication is that $`R_{\mathrm{in}}`$ increases by a factor of $``$110. The power spectra we measure for XTE J1650–500 at low luminosities are consistent with the standard broken power-law or zero-centered Lorentzian models that are typically used to fit Hard state power spectra (Tomsick et al., 2003b), and it is interesting that optical observations of other black hole systems in quiescence show similar power spectra but with even lower break frequencies (see Hynes et al., 2003, and the paper by Shahbaz et al. in these proceedings). Although we are comparing optical and X-ray properties, the similar timing properties suggest that black hole transients have the same overall structure and that the same physical processes are operating in the Hard state and in quiescence (Tomsick et al., 2003b). Furthermore, the fact that the quiescent energy spectra of black holes exhibit a component at optical and UV energies that may be caused by thermal emission from a truncated disk (McClintock et al., 2003) underscores the idea that physical properties change gradually as the source luminosity drops in the Hard state and quiescence but that the overall structure does not undergo a major transition.

## 6 Theoretical Models for the Hard State

Table 1 compares the predictions of three different theoretical pictures for the Hard State to observational properties. In the “Sphere+Disk” model, the accretion disk is truncated at some value of $`R_{\mathrm{in}}`$ and a spherical corona is present for radii closer to the black hole (Dove et al., 1997; Esin et al., 1997). The “Slab” model is typically envisioned as an optically thick disk extending close to the ISCO with active regions above the disk acting as the corona to produce the hard X-ray emission (Galeev et al., 1979). While the emission mechanism for these first two models is inverse Componization, the third, which is the X-ray jet model, produces emission via a synchrotron mechanism (Markoff et al., 2001).

With the additional free parameter, $`R_{\mathrm{in}}`$, the Sphere+Disk (SD) model provides a natural explanation for the reflection strength vs. $`\mathrm{\Gamma }`$ and characteristic frequency vs. $`\mathrm{\Gamma }`$ correlations described above. In addition, the luminosity independence of states argues for the addition of a parameter that may not be directly tied to mass accretion rate. However, while there is observational evidence that jets are present in the Hard state222Steady, compact radio jets have been resolved in the Hard state for GRS 1915+105 and Cyg X-1 (Dhawan et al., 2000; Stirling et al., 2001). The radio emission from other Hard state systems is likely also from unresolved jets., jets have not been incorporated into SD models such as the Advection-Dominated Accretion Flow (ADAF) model. It has previously been suggested that accretion energy in the Hard state might be used to power a jet (Blandford and Begelman, 1999; Fender et al., 2003), but it is not clear if, e.g., the Advection-Dominated Inflow-Outflow Solution (ADIOS) can explain the Hard state radio emission.

For the Slab model, it has been demonstrated in more detail that a magnetic corona can drive a jet (Merloni and Fabian, 2002). However, in this model, explanations for the reflection and frequency vs. $`\mathrm{\Gamma }`$ correlations are unclear. Perhaps the size of the corona changes, but more work is necessary to determine if the correlations can be reproduced within this model. A significant point in favor of the X-ray jet model is that it can explain, in detail, the reported correlations between X-ray and radio flux (see Gallo et al., 2003, and the paper by Corbel in these proceedings); however, similar X-ray/radio flux correlations may arise in more general situations with non-X-ray-emitting jets (Heinz and Sunyaev, 2003). More work is necessary to determine if a synchrotron X-ray jet is viable. As shown in Zdziarski et al. (2003), a synchrotron cutoff is expected to be more gradual than the exponential cutoffs that are typically observed in black hole spectra. Also, as discussed in Tomsick et al. (2003b), it is not clear if the X-ray jet model can reproduce a change in the characteristic frequencies of the system by a factor of 1200 since this would seem to imply a very large change in the location of the shock that produces the high energy electrons that give rise to X-ray emission.

## 7 Beyond RXTE

Studies of black hole X-ray spectra require a broad bandpass as well as flexible scheduling due to the rapid evolution that these sources can exhibit. With the combination of an all-sky monitor, instrumentation (PCA and HEXTE) covering the 3-200 keV energy range, and a dedicated Project Scientist and scheduling team, RXTE has been an excellent tool for research on black hole spectra.

An instrument with an even larger effective area, as has been discussed at this meeting, would be useful for spectral studies. Black hole systems are highly variable, and a larger area would allow for high quality spectra to be obtained over a shorter period of time. For example, it may be possible to study spectral evolution on sub-second time scales when QPOs are present. Also, a larger effective area above $``$50 keV would allow for detailed studies of the evolution of the high energy component. Although I have not discussed measurements of time lags and coherence, our group’s observations of sources during state transitions indicate that these parameters provide useful information on the physics of the system (Kalemci et al., 2003), and a larger area would be very beneficial to such measurements.

The spectra shown in Figure 2 show that the RXTE bandpass captures the two primary spectral components, but it should be noted that important information would be lost with a reduced bandpass. Background is another consideration as the transition to the Hard state during outburst decay can occur at 3-25 keV flux levels below $`1\times 10^9`$ erg cm-2 s-1 even for relatively nearby systems with distances of a few kpc. A significantly higher background level could limit black hole studies to the brightest portions of the outburst. Finally, as most X-ray binaries are in the Galactic plane, source confusion is a concern for studying fainter sources, and this can be problematic for RXTE, given the $`1^{}`$ radius (FWZI) field of view. A smaller field of view could be beneficial for studies of fainter sources in the Galactic plane.

I would like to thank my RXTE collaborators Emrah Kalemci, Stephane Corbel, and Philip Kaaret, and I am especially grateful to Jean Swank for all of her hard work and useful suggestions while planning and co-ordinating our rather demanding observing program. I acknowledge partial support from NASA grant NAG5-13055.


format: text

text:
# Adsorption of Xe and Ar on Quasicrystalline Al-Ni-Co

## I Introduction

The growth and equilibrium structure of an adsorbed film are dependent on the competing adsorbate-adsorbate and adsorbate-substrate interactions. The laterally aperiodic adsorption potential of a quasicrystalline (QC) surface provides an interesting case of competing interactions for rare gas adsorbates, which favor close-packed monolayer structures in the absence of substrate corrugation. Such incongruity often produces new phenomena which are interesting and exotic on a fundamental level. In this case, because quasicrystals can have radically different physical properties than their periodic counterparts, it also provides a significant tool for the design and growth of thin films having specific properties.

The 10-fold surface of decagonal Al-Ni-Co is aperiodic in the surface plane but periodic in the perpendicular direction Steurer ; cervellino ; takakura ; saitoh ; steinhardt ; haibach ; mihalk ; henley ; cockayne ; abe ; zaharko . Each quasicrystalline layer comprises an aperiodic array of Al, Ni and Co with points of 5-fold rotational symmetry. Each plane is related to its neighboring planes by a rotation of 36 degrees, producing an ABAB stacking sequence. The structure of this surface has been studied using various techniques, including low-energy electron diffraction (LEED) and scanning tunneling microscopy (STM). A combination of these techniques was used recently to demonstrate that the surface structure of Al-Ni-Co is similar to the bulk structure determined by x-ray diffraction Steurer , but with a significant relaxation of the top layer and some intralayer buckling Nicola .

There have been many theoretical studies of physical adsorption on flat substrates and on periodic substrates Milton ; landoldt , but none on aperiodic substrates. Recently a rigid-lattice total energy calculation for Al adsorption on a quasiperiodic substrate produced some interesting results concerning the growth and size of Al clusters on the surface lusher . That calculation was carried out by assuming a Lennard-Jones (LJ) potential between the adsorbate Al atoms and the substrate atoms, assuming a bulk-like structure for the surface. Comparisons were made to the results of LEED experiments for the structure, orientation and domain size of the Al film. In particular, the domain size distribution was shown to be a direct consequence of the competing interactions in the substrate and the film.

Because the interactions are weaker, simpler and better known for physisorbed gases, we believe we can gain a fuller understanding of the effects of aperiodicity and symmetry by studying rare gas adsorption on Al-Ni-Co. In this study, we have chosen Xe and Ar as the adsorbates. We have calculated the gas-surface adsorption potential and the adsorption properties in the low-coverage limit using a virial expansion. The results of these calculations are compared to the results for adsorption on periodic substrates, and to thermodynamic LEED measurements of Xe adsorption on quasicrystalline Al-Ni-Co.

## II Experimental Characterization of the Al-Ni-Co Substrate and Xe Adsorption

The decagonal Al-Ni-Co quasicrystal has a basic structure that consists of a stack of identical or nearly-identical 5-fold symmetric planes, each related to adjacent ones by a $`\pi /5`$ rotation cervellino . This produces a stacking structure (ABAB) of A and B terminations, with a 10-fold screw axis. A schematic representation of the atomic positions of one of such planes is shown in Figure 1, where the different chemical identities and different local geometry among the same chemical identity are specified. The 10-fold surface presents a contraction (10$`\%`$) of the first layer and an expansion (5$`\%`$) of the second layer, with a small degree of intralayer rumpling; in plane reconstruction is minimal, if present. The two-dimensional atomic density of a layer is 0.123 $`\stackrel{̊}{A}`$ nicola2 .

The quasicrystalline d-$`Al_{73}Ni_{10}Co_{17}`$(Co-rich phase) sample was grown at Ames National laboratories using the decantation method moog . The surface was obtained by a cut perpendicular to the 10-fold axis, and then polished as described elsewhere nicola2  to obtain a surface that was within 0.5  of the ten-fold orientation. The sample preparation in ultra-high vacuum consisted of several cycles of Ar ion bombardment (500 eV ions) for about 45 minutes, followed by annealing for six to eight hours at temperatures up to 1060 K, as measured by a K-type thermocouple in contact with the sample, and an optical pyrometer. The LEED pattern after preparation was observed to have well-defined spots and relatively low intensity between the primary spots, as shown in Figure 2. The impurity level was below detectablility, as measured by Auger electron spectroscopy. The symmetry of the LEED pattern is 10-fold, due to the presence of two equivalent surface terminations rotated by 36  renee1 . This method of preparation has been shown to produce a surface having a structure that is essentially identical to the bulk structure with the exception of a small degree of surface relaxation Nicola . The LEED intensities of the diffraction spots were measured using a rear-view LEED system, with the electron beam at normal incidence to the surface. The LEED adsorption isobars were obtained by holding the Xe pressure at a fixed value while changing the temperature and acquiring LEED frames. For each frame the integrated spot intensity of a substrate diffraction peak was extracted and plotted versus the acquisition temperature, giving the adsorption isobar nicola2 . The temperature was measured using a chromel-alumel thermocouple in contact with the sample, in a range between 60 and 140K. The coverage is assumed to be linearly related to the attenuation of the intensity of non-specular substrate diffraction peak in the submonolayer range, where one monolayer is defined by the break in the isobar. Thus, the substrate peak intensity before adsorption is defined to correspond to zero coverage, the intensity at the break in the isobar is defined to be one monolayer, and the intermediate coverages vary inversely and linearly with the intensity moog .

## III Calculation of the QC interaction potential

The calculation of the interaction energy between Xe (or Ar) and the Al-Ni-Co quasicrystalline structure is performed by summing over the LJ pair interactions, Xe-Al, Xe-Ni and Xe-Co. In the following, we will refer to Xe as the adsorbent, but a similar technique is used for Ar for which we will also present results. Lusher and coworkers employed a similar LJ superposition while exploring the growth of Al on Al-Co-Ni QC surface lusher . They assumed that the interaction of Al with each species of the QC is the same. Using this approximation, they were able to predict nanocrystal growth of Al on the QC substrate. Here, we assume that there are two distinct pair interactions: Xe-Al and Xe-TM, where TM stands for either transition metal, Co or Ni.

In order to find the Xe-Al and Xe-TM LJ parameters, we consider adsorption of Xe on separate elemental surfaces of Al(110) and Ni(100). The semiempirical arithmetic combining rules between a Xe gas ($`\sigma _{Xe}=4.1\stackrel{̊}{A}`$) and Al or Ni crystals ($`\sigma _{Al}=2.5\stackrel{̊}{A},\sigma _{Ni}=2.2\stackrel{̊}{A}`$) yield $`\sigma _{XeAl}=3.3\stackrel{̊}{A}`$ and $`\sigma _{XeTM}=3.1\stackrel{̊}{A}`$. In order to estimate the gas-solid well-depth $`ϵ`$ parameters, we use the experimental heats of adsorption of Xe on Al(110) ($`Q_{XeAl}=190`$ meV) and on Ni(100) ($`Q_{XeNi}=226`$ meV) landoldt , which are basically a measure of the interaction well depth, and the summed pairwise interaction of Xe with Al(110) and Ni(100) crystals, respectively steele . Finally, we find $`ϵ_{XeAl}24`$ meV and $`ϵ_{XeTM}23`$ meV. Notice that the LJ parameters of Xe-Al and Xe-TM are not very different. However, the small difference in the parameter values and the small rumpling of the QC surface play a role in the adsorption potential, as we will show in the following.

Once the Xe-Al and Xe-TM interactions are obtained, we sum them to produce the net Xe-QC interaction. Before discussing the results, we note that they should be interpreted with caution since summation of pair potentials omits effects of electron delocalization, which have been found to occur on metallic substrates renee .

As described in the introduction, the coordinates of the Al-Ni-Co QC atoms were found using LEED data Nicola . Due to the QC’s aperiodicity, periodic boundary conditions are not appropriate and all substrate atoms have to be considered in calculations. We simplify by selecting a computational cell that is large enough ($`56\stackrel{̊}{A}\times 56\stackrel{̊}{A}`$) to describe the potential adequately. One way to depict the potential energy $`V(𝐫)`$ is to construct a “minimum energy surface”, defined as follows: for any surface-parallel position $`𝐑`$ =(x,y) we evaluate the position z$`(𝐑)`$=$`z_{min}(𝐑)`$ at which $`V(𝐫)`$ is a minimum. The resulting well-depth and force constant are called $`D(𝐑)`$ and $`k(𝐑)`$, respectively.

A contour plot of the minimum energy surface is shown in Figure 3(a), where the various curves represent isopotential surfaces. The minimum energy surface exhibits the five-fold symmetry of the QC top layer. The big (red) circles correspond to the repulsive regions of the minimum energy surface and are placed on top of the Al-2 species (see Fig. 1). In general, the repulsive part of the minimum energy surface follows closely the distribution of atoms of the top layer, with one exception: the TM atoms which are buried deepest in the first layer. The insensitivity of the interaction energy to those atoms is due to the fact that $`\sigma _{XeTM}`$ and $`ϵ_{XeTM}`$ are smaller than $`\sigma _{XeAl}`$ and $`ϵ_{XeAl}`$.

The most striking feature of the minimum energy surface is its high “corrugation”: local potential minima range from -150 meV to -270 meV. A corrugation of this order was also found on open crystalline structures such as reconstructed Si and Ge conrad ; packard , but is not usual for metallic surfaces. However, this corrugation does not have the meaning of the hopping barrier between adjacent sites. To better illustrate this point, we plot the lateral variation of $`V_{min}(x,y)`$ along different path in Figure 3 (b). The full and dashed curves correspond to vertical paths in Fig. 3(a) at x=-1.6 $`\stackrel{̊}{A}`$ and x=1.6 $`\stackrel{̊}{A}`$. The energy difference between local minima and bridges is very irregular and ranges in general from few meV to $``$ 100 meV. We attribute this large range of local minima to the QC’s aperiodicity as well as the “holes” ( 4 $`\stackrel{̊}{A}`$ width) of its lateral structure Nicola . A more complete analysis of the hopping barriers between adjacent sites, with applications in diffusion, involves an identification of the adsorption sites and the saddle points between those. This investigation is in progress.

In order to illustrate the most attractive adsorption sites, in Figure 4 we plot the regions where the adsorption potential has values between -270 meV and -220 meV. McGrath et al argued that quasicrystals could be potentially used as templates for quasicrystalline 2D structures mcgrath1 ; mcgrath2 . With this in mind, they have studied adsorption of $`C_{60}`$ on the Al-Pd-Mn quasicrystal surface, which is known to exhibit depressions of about 7 $`\stackrel{̊}{A}`$ width in a pentagrid structure. STM images at low coverage showed that $`C_{60}`$ molecules occupy some of those pentagonal holes. In a similar way, Al-Ni-Co is a good candidate as a template for Xe or Ar. It is not clear from Figure 4 if adsorption of Xe at low coverage would lead to a quasicrystalline 2D long-range order structure. Monte Carlo simulations are in progress to investigate this. The interaction energy between Ar and QC shows very similar properties, the main difference being that the local potential minima range from -60 meV to -130 meV.

Another way to depict this highly corrugated potential is provided in Figure 5(a). There is shown the lateral variation of $`V(𝐫)`$ in the x-z plane for constant z and y=0 (see Figure 3 a). As expected, far from the surface (z=5.8 $`\stackrel{̊}{A}`$), the potential is high and nearly independent of x. Close to the surface (z=2.8 $`\stackrel{̊}{A}`$), the potential exhibits a large and irregular corrugation. An unexpected feature is the following: at some particular values of x (x $``$ -25, - 15, -10, -2.5, 5, 20 $`\stackrel{̊}{A}`$) the potential is “more” attractive close to the surface (z = 2.8 $`\stackrel{̊}{A}`$), but is “less” attractive at distances farther away from the surface (z = 3.8, 4.4 $`\stackrel{̊}{A}`$). Recently, it has been argued that the corrugation can change sign even for an elemental crystal jean . This may lead to a possible steering effect on adsorption on the QC, especially important in the process of film growing.

To better understand the interaction potential landscape, we define a “volume density of states” f(V) as follows: f(V)dV equals the volume above the QC surface such that the potential energy lies in the interval \[V,V+dV\]. Then

$$f(V)=𝑑𝐫\delta [VV(𝐫)]=𝑑𝐑\left\{\frac{1}{|V/z|^{(1)}}+\frac{1}{|V/z|^{(2)}}\right\}\theta [V+D(𝐑)]$$
(1)
where $`\theta (x)`$ is the unit step function. In the second expression, the denominators equal the magnitudes of the surface-normal forces $`|V/z|^{(1)}`$ and $`|V/z|^{(2)}`$ evaluated at the two points (1,2) where $`V=V(𝐫)`$. Equation 1 has been numerically integrated and the result appears as the full curve in Figure 5(b). Note that the full curve exhibits a broad maximum around -150 meV, due to the heterogeneity of the substrate. We may understand this by considering the behavior of f(V) in two limits, low V and high V. The region of the potential energy minimum (low V) is particularly important. At a given $`(𝐑)`$, this may be evaluated with a local harmonic approximation : $`V(𝐫)D(𝐑)+k(𝐑)[zz_{min}(𝐑)]^2/2`$ and $`|V/z|^{(1)}|V/z|^{(2)}=[2k(𝐑)(V+D(𝐑))]^{1/2}`$. Therefore, f(V) becomes:

$$f_{har}(V)=𝑑𝐑\theta [V+D(𝐑)]\sqrt{\frac{2}{k(𝐑)(V+D(𝐑))}}$$
(2)
The local harmonic approximation for the heterogeneous surface is shown as a dashed curve in Figure 5(b). The agreement with the true f(V) is excellent at low energy, reproducing the broad maximum due to the heterogeneity.

The harmonic approximation to f(V) can be further approximated by assuming that the behavior in the region of the potential energy minimum is the same for all $`𝐑`$, $`D(𝐑)=D`$ and $`k(𝐑)=k`$. This corresponds to a monolayer film on a flat surface, for which we choose as well depth the average value $`D_{ave}(𝐑)`$=180 meV. In this case:

$$f_{flat}(V)=A\theta (V+D)\sqrt{\frac{2}{k(V+D)}}.$$
(3)
Here A is the computational surface area. Results of Equation 3 are shown in Figure 5(b) as the dotted curve. Notice that instead of a shoulder, $`f_{flat}(V)`$ exhibits a singularity at $`V=D_{ave}`$; manifestly, the smooth surface approximation is not suitable for a QC.

At high energy, another approximation can be employed, which is related to the attractive part of the LJ potential. For a semi-infinte substrate, at positions far from the surface $`V(z)=C_3/z^3`$, where $`C_3=(\pi /6)_{i=1}^\gamma n_iC_6^{(i)}`$ for a substrate composed of $`\gamma `$ species of atoms, and $`C_6^{(i)}=4ϵ_i\sigma _i^6`$ is the coefficient of the attractive part in the usual LJ potential. Using this potential in Equation 1, we obtain:

$$f_{high}(V)=\frac{A}{3}(\frac{C_3}{|V|^4})^{1/3}$$
(4)
Equation 4 is shown as dot-dashed curve in Figure 5(b). The approximation employed at high energy (dot-dash) does not match exactly the true f(V) (full curve) since it assumes a semi-infinite substrate (which is not the case of the QC used in calculations) and since the volume considered in the numerical calculation extends to $`z=\mathrm{}`$ while in practice the numerical calculation stopped at $`z_{max}=10\stackrel{̊}{A}`$.

## IV Adsorption at low coverage

Adsorption at low coverages can be described in terms of a virial expansion steele ; Milton .

$$N_{excess}=k_H\beta P+C_{AAS}\beta ^2P^2+\mathrm{}$$
(5)
The leading term of this expansion is all that is needed when the coverage is so low that interactions between adsorbate atoms can be neglected. In this regime, called the Henry’s law regime, the coverage is proportional to the pressure (P) and the coefficient of proportionality is called the Henry’s law constant:

$$k_H=d^3𝐫[e^{\beta V(𝐫)}1].$$
(6)
Here $`\beta =1/(k_BT)`$ and $`V(𝐫)`$ is the adsorbate-substrate interaction potential. At somewhat higher coverage, the interactions between adatoms become important and are taken into account through the second virial coefficient $`C_{AAS}`$:

$$C_{AAS}=d^3𝐫_\mathrm{𝟏}d^3𝐫_\mathrm{𝟐}e^{\beta [V(𝐫_\mathrm{𝟏})+V(𝐫_\mathrm{𝟐})]}[e^{\beta u(|𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐}|)}1]$$
(7)
where $`u(|𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐}|)`$ is the interaction between two adatoms. While the Henry’s law coefficient involves a 3D integral which can be relatively easily calculated, the second virial coefficient is a 6D integral which is computationally expensive due to the large domain of integration.

Figure 6(a) presents an Arrhenius plot of the coverage (or excess number per unit area) as a function of 1/T for Xe and Ar in two regimes: very low coverage (Henry’s law regime) and higher coverage ($`C_{AAS}`$ included). At high T, the contribution of the second virial term is negligible, whereas at low T, mutual adatom interactions become important and the coverage is enhanced by the $`C_{AAS}`$ term. However, when the second virial term greatly exceeds the first term, the virial expansion becomes divergent and its results cannot be trusted. Figure 6(b) presents a comparison of our calculated results with an experimental isobar. The agreement between calculation and experiment is rather good at low coverage (below $`0.01A^2`$) considering the simplicity of the model.

An interesting question arises in the heterogeneous environment of the QC potential: how does this heterogeneity affect the mutual adatoms’ interactions? To investigate this, we compute the second virial coefficient in the 2D approximation steele :

$$B_{2D}=\frac{C_{AAS}}{2k_H^2}$$
(8)
On a substrate modelled as a flat 2D continuum, a monolayer film is perfectly mobile, so $`B_{2D}`$ is not influenced by substrate but only by the adatoms’ interaction $`u(r)`$:

$$B_{2D}^{flat}=1/2𝑑r[e^{\beta u(r)}1]$$
(9)
Figure 7 presents a comparison between $`B_{2D}`$/area for a flat substrate and for the QC. Note that the QC’s corrugation enhances the effect of the adatoms’ attraction since $`B_{2D}`$ is larger in magnitude for a QC than for a flat substrate. We believe that this is another consequence of the aperiodicity and the semi-close packed structure of the QC. The increased attraction (also evident in the isosteric heat) occurs because the corrugated adsorption potential tends to bring particles closer together than a smooth surface. This effect is much stronger than in the case of a crystal, where the regular corrugated potential enhances the interactions’ effect to a smaller extent steele .

The isosteric heat of adsorption can be computed from the equation of state:

$$Q_{st}=(\frac{d\mathrm{ln}P}{d\beta })_N$$
(10)
In the Henry’s law regime, the isosteric heat reduces to:

$$Q_{st}=k_BT<V>$$
(11)
where $`<V>=(𝑑rVe^{\beta V})/(𝑑re^{\beta V})`$ is the mean value of the adsorption potential. In the case of a flat substrate, assuming a harmonic adsorption potential, $`<V>=D+k_BT/2`$. Thus, the isosteric heat of a perfectly mobile monolayer on a flat surface increases monotonically with T as $`Q_{st}=D+k_BT/2`$. For a heterogeneous substrate, however, the dependence of $`Q_{st}`$ on T must be determined numerically. Figure 8 presents the isosteric heat results for Xe and Ar in the Henry’s law regime (coverage independent) and including the second virial coefficient correction (coverage dependent). Note that in both cases $`Q_{st}Dk_BT/2`$, meaning that $`<V>D+3k_BT/2`$; this occurs because Xe atoms tend to be confined three-dimensionally on the QC (rather than 1D in the smooth surface mode). This finding is consistent with the very corrugated potential in Figure 3(b). A comparison between the Henry’s law regime and second virial coefficient correction shows that the effect of interactions is to increase the isosteric heat. At T= 70 K, the calculated isosteric heat of Xe is 290 meV in the Henry’s law regime and 305 meV when interactions are included. Both predictions are close to the value of 300 meV at 0.25 monolayer adsorption extracted from experimental isobarsNicola .

## V Conclusions

In conclusion, we have computed the adsorption potential of Xe and Ar adsorbates on Al-Ni-Co quasicrystalline substrate by adding gas-surface LJ interactions. Due to the substrate’s aperiodicity and heterogeneity, the adsorption potential was found to be highly corrugated, with local minima ranging from -150 to -270 meV for Xe and -60 to -130 meV for Ar. The minimum energy surface of the adsorption potential exhibits the 5-fold symmetry of the QC’s top layer, with the most repulsive regions located on top of the QC’s atoms and the most attractive regions in the “holes” of the semi-close packed structure of the QC.

To understand the distribution of adsorption potential values, we have defined a volume density of states f(V) and computed this quantity for QC, for a flat surface and for a harmonic approximation of the interaction potential. While $`f_{flat}(V)`$ exhibits a singularity when $`V=D`$ (the well depth), $`f_{QC}(V)`$ has a broad maximum due to the QC’s “heterogeneity”. This shoulder is well reproduced by a local harmonic approximation of the adsorption potential in each point above the QC surface.

To model adsorption at low coverages, a virial expansion of the equation of state was employed, including first and second virial coefficient terms. Isobars obtained from the virial expansion agree well with experimental isobars of Xe on Al-Ni-Co QC. The isosteric heat of Xe was found to be 305 meV, close to experimental value of 300 meV at 0.25 monolayer. The dependence of the isosteric heat on T indicates that adatoms are 3D-confined in the attractive regions of the adsorption potential. This finding is consistent with the results of the second virial coefficient which suggest that the corrugation of the potential tends to bring adatoms together and enhance their mutual interactions. Future work will focus on a Monte Carlo study at arbitrarily high coverage and on a search for the adsorbate ground state structure.

## VI Acknowledgments

We thank L.W. Bruch and G.D. Mahan for edifying comments and the National Science Foundation (Grants 02-08520 and 03-03916) for support of this research.

## VII Figure captions

1. One layer of a $`45\times 45\stackrel{̊}{A}^2`$ slab, representing one 5-fold plane of the d-AlNiCo quasicrystal surface, along the 10-fold axis. The atom colors correspond to sublayer groups differentiated according to the chemical identity and the local geometry (nearest neighbor distance): Red - Transition Metals (Ni or Co) 1, Black - Transition Metals 2, Green Al-1, Blue Al-2, Cyan Al-3, Yellow Al-4 Nicola .

2. LEED pattern from the 10-fold surface of the clean d-AlNiCo, acquired at 60 K. The incident energy was 55 eV.

3. Interaction potential energy between Xe (or Ar) and Al-Ni-Co QC: (a) Contour plot of the minimum energy surface. Different colors correspond to different isopotential curves: red = -175 meV, green = -200 meV, turqoise = -225 meV, blue = -250 meV; (b) Lateral variation of the minimum potential along two paths: x= 1.6 $`\stackrel{̊}{A}`$ (dashed curve) and x=-1.6 $`\stackrel{̊}{A}`$ (full curve).

4. The preferred adsorption regions at low coverage (-270 meV $`V_{min}(r)`$ -220 meV).

5. Interaction potential energy between Xe and Al-Ni-Co QC: (a) Potential energy as a function of x for constant z above the QC, (b) The “volume density of states” of the QC (full curve), harmonic approximation limit at low energy (dashed), the high energy limit approximation (dot-dashed) and for a flat surface (dotted).

6. Coverage as a function of temperature (a) Arrhenius plot of Xe and Ar isobars at P=1.6\*$`10^7`$ mbar in the Henry’s law regime (Ar-full, Xe-dot) and including the second virial term (Ar-dot-dash, Xe-dash), (b) Comparison of calculations with an experimental isobar of Xe at low coverage.

7. The second virial coefficient for a flat surface and the QC surface.

8. Isosteric heat of Xe and Ar on the QC in the Henry’s law regime: Xe (dashed), Ar (dotted), and including the second virial coefficient correction: Xe (full), Ar (dot-dashed).


text:
# Lattice gas cellular automata model for rippling and aggregation in myxobacteria

## 1 Introduction

Myxobacteria are one of the prime model systems for studying cell-cell interaction and cell organization preceding differentiation. Myxobacteria are social bacteria which swarm, feed and develop cooperatively . When starved, myxobacteria self-organize into a three-dimensional fruiting body structure. Fruiting body formation is a complex multi-step process of alignment, rippling, streaming and aggregation that culminates in the differentiation of highly elongated, motile cells into round, non-motile spores. A successful model exists for the fruiting body formation of the eukaryotic slime mold Dictyostelium discoideum . Understanding the formation of fruiting bodies in myxobacteria, however, would provide a new insight since collective myxobacteria motion depends not on chemotaxis as in Dictyostelium but on contact-mediated signaling (see  for a review).

During fruiting body formation myxobacteria cells are elongated, with a 10:1 length to width ratio, and move along surfaces by gliding. Gliding occurs in the direction of a cell’s long axis  and is controlled by two distinct motility systems in myxobacteria . One of the most interesting patterns that develops during myxobacteria morphogenesis is rippling, which often occurs spontaneously and transiently during the aggregation phase . Rippling myxobacteria form equidistant ridges of high cell density which appear to advance through the population as rhythmically traveling waves  (Figure 1). Cell movement in a ripple is approximately one-dimensional since the majority of cells are aligned and move in parallel lines with or against the direction of wave propagation . Tracking individual bacteria within a ripple has shown that cells reverse their traveling directions back and forth and that each travels on the order of one wavelength between reversals . The ripple waves propagate with no net transport of cells  and wave overlap causes neither constructive nor destructive interference . Although mechanisms for gliding are not fully understood, they are believed to account for both alignment  and reversals  in myxobacteria.

Rippling is related to a membrane-associated signaling protein called C-factor. C-factor regulates rippling , cells without the ability to produce C-factor fail to ripple  and the addition of C-factor (extracted from fruiting body cultures) causes cell reversal frequencies to increase three-fold . C-signaling occurs via the direct cell-cell transfer of C-factor when two elongated cells collide head to head . Understanding the mechanisms of the rippling phase may reveal many clues about the way myxobacteria organize collective motion since C-factor is also involved in all other stages of fruiting body formation. For example, cells lacking in C-factor fail to aggregate or sporulate  while high concentrations of exogeneous C-factor induce aggregation and sporulation .

In this paper we use two lattice gas cellular automata (LGCA) models to simulate rippling and aggregation during the fruiting body formation of myxobacteria, to show the potential of cellular automata as models for biological pattern formation processes, and to evaluate, in particular, the necessity of different biological assumptions shown in previous models for pattern formation in myxobacteria.

Sager and Kaiser  have proposed that precise reflection explains the lack of interference between wave-fronts for myxobacteria rippling. Oriented collisions between cells initiate C-signaling that causes cell reversals. According to this hypothesis of precise reflection, when two wave-fronts collide, the cells reflect one another, pair by pair, in a precise way that preserves the wave structure in mirror image. Figure 2 shows a schematic diagram of this reflection.

We present a new LGCA approach for modeling cells which is computationally efficient yet approximates continuum dynamics more closely than assuming point-like cells. As an example of this new approach, we present a model for myxobacteria rippling based on the hypothesis of precise reflection and a model for aggregation based on C-signaling.

This paper is organized as follows. The biological assumptions for precise reflection and C-signaling that motivate the models are described in the next section. In section 3 we describe specifics of two LGCA models. In section 4 results of modeling rippling phenomenon are discussed in detail. Section 5 provides description of a model for aggregation centers. The paper ends with a summary section.

## 2 Biological Background

In this section we describe the biological observations which motivate our models for rippling and aggregation.

Rippling and aggregation are both controlled by C-signaling and are characterized by specific high cell density patterns (in particular, moving high density ridges in rippling and stationary high density mounds in aggregation). There is a marked relationship between cell density, levels of C-signaling and behaviors in myxobacteria triggered by C-signaling . C-signaling increases with density since end-to-end contacts between cells are more likely with increased density  and high cell densities favor spatial arrangements in which there are many end-to-end contacts due to the polarity of myxobacterial cells . Cell density and C-signaling levels increase together from rippling to aggregation and from aggregation to sporulation . Further, increased thresholds of C-factor induce rippling, aggregation and sporulation respectively , suggesting C-signaling levels, as a measure of cell density, are checkpoints for different stages of development. Kim and Kaiser suggest that C-factor may act as a developmental timer that triggers sporulation only when cell density is as high as possible . A high density aggregate will culminate in a fruiting body with a large number of spores ensuring that the next cycle is started by a population of cells . Sager and Kaiser have also observed the effect of C-signaling-competent cell density upon ripple wavelength . They dilute a cell population of C-signaling-competent cells with cells that are able to respond to C-factor but are not able to transmit it. They find that with increased concentrations of these csgA-minus cells, ripple wavelength increases non-linearly.

In addition to cell density patterns, cell geometries are important throughout the stages of fruiting body formation and distinguish different stages. During the fruiting body formation, cells form aligned patches from a random distribution . For rippling, a large number of cells must be aligned both parallel and anti-parallel within the same field. For streaming, cells form long chains which flow cooperatively in aggregation centers . In Stigmatella spp., cells moving in circles or spirals form microscopic transient aggregates. These aggregates disappear as cells also spiral away tangentially . Macroscopic aggregates form in areas of high density  and may also disappear as cells apparently stream along chains from one aggregate to the other . The mature structures of fruiting bodies are diverse and species-dependent, ranging in size between 10 and 1000 $`\mu `$m . In Myxococcus xanthus, the basal region of the fruiting body is a shell of densely packed cells which orbit in two directions, both clock-wise and counter-clockwise, around an inner region only one-third as dense . In Stigmatella aggregates, cells are organized in concentric circles or ellipses and cells move in a spiral fashion up the aggregate as the fruiting body develops .

Current models for rippling () assume precise reflection. Key differences among these models include their biological assumptions regarding the existence of internal biochemical cell cycles. It is still not known if an internal cell timer is involved in myxobacterial rippling. Several models with completely different assumptions all qualitatively produce ripple patterns resembling experiment.

An internal timer is a hypothetical molecular cell clock which regulates the interval between reversals. The internal timer may specify a delay, or minimum period between reversals, which would include the refractory period, see below, and a minimum response time; the minimum period of time required for a non-refractory cell to become stimulated to turn. Also, the internal timer may specify a maximum oscillation period, in which case the timer may speed up or slow down depending upon collisions, but the cell will always turn within a specified period of time even without collisions. Individual pre-rippling cells reverse spontaneously every 5-10 minutes with a variance in the period much smaller than the mean . This would suggest that there is a component of the timer specifying a maximum oscillation period. Also, observation of rippling bacteria reveals that cells oscillate even in ripple troughs where the density is too low for frequent collisions  further supporting the hypothesis of a maximum oscillation period.

The refractory period is a period of time immediately following a cell reversal, during which the cell is insensitive to C-factor. Although there is no evidence of a refractory period in the C-signaling system, the refractory period is a general feature of bacterial signaling systems  (for a description of the role the refractory period plays in Dictyostelium, see ). The addition of .02 units of external C-factor triples the reversal frequency of single cells from .09 reversals per minute to .32 reversals per minute . Cells do not reverse more frequently at still higher levels of C-factor, however, suggesting the existence of a minimum oscillation period of 3 minutes in response to C-factor. This minimum oscillation period would be the sum of the refractory period and the minimum response, so the duration of the refractory period cannot be guessed from this fact alone.

To resolve the conflicts of these models for rippling our first LGCA model is designed to test different assumptions. The results of our model for rippling shows that rippling is stable for a wide range of parameters, C-signaling plays an important role in modulating cell density during rippling, and non-C-signaling cells have no effect on the rippling pattern when mixed with wild-type cells. Further, by comparing model results with experiments, we can conclude reversals during rippling would not be regulated by a built-in maximum oscillation period.

We then present a second LGCA model for aggregation based on C-signal alignment, which reproduces the sequence and geometry of the non-rippling stages of fruiting body formation in detail, showing that a simple local rule based on C-signaling can account for many experimental observations.

## 3 Model and Method

LGCA are relatively simple Cellular Automata models. They employ a regular, finite lattice and include a finite set of particle states, an interaction neighborhood and local rules that determine the particles’ movements and transitions between states . LGCA differ from traditional CA by assuming particle motion and an exclusion principle. The connectivity of the lattice fixes the number of allowed non-zero velocities or channels for each particle. For example, a nearest-neighbor square lattice has four non-zero allowed channels. The channel specifies the direction and magnitude of movement, which may include zero velocity (resting). In a simple exclusion rule, only one particle may have each allowed non-zero velocity at each lattice site. Thus, a set of Boolean variables describes the occupation of each allowed particle state: occupied (1) or empty (0). Each lattice site on a square lattice can then contain from zero to four particles with non-zero velocity.

The transition rule of an LGCA has two steps. An interaction step updates the state of each particle at each lattice site. Particles may change velocity state, appear or disappear in any number of ways as long as they do not violate the exclusion principle. In the transport step, cells move synchronously in the direction and by the distance specified by their velocity state. Synchronous transport prevents particle collisions which would violate the exclusion principle (other models define a collision resolution algorithm). LGCA models are specially constructed to allow parallel synchronous movement and updating of a large number of particles .

### 3.1 Representation of cells

In classical LGCA, biological cells are dimensionless and represented as a single occupied node on a lattice (e.g., see  and ). Interaction neighborhoods are typically nearest-neighbor or next-nearest-neighbor on a square lattice. The exclusion principle makes transport unwieldy when a single cell occupies more than one node since a cell may only advance if all the channels it would occupy are available. Similarly, it is difficult to model the overlapping and stacking of cells. Cells without dimension are untenable for a sophisticated model of myxobacteria fruiting body formation, however. Cells are very elongated during rippling, streaming and aggregation and form regular, dense arrays by cell alignment. Also, a realistic model of cell overlap and cell stacking is needed since interaction occurs only at specific regions of highly elongated cells and cell density is a critical parameter throughout this morphogenesis.

Börner et al.  have mediated the problem of stacking by introducing a semi-three-dimensional lattice where a third z-coordinate gives the vertical position of each cell when it is stacked upon other cells. Stevens  has introduced a model of rod-shaped cells that occupy many nodes and have variable shape in her cellular automata model of streaming and aggregation in myxobacteria. Neither of these two models are LGCA since they do not incorporate synchronous transport along channels. We device a novel way of representing cells which facilitates variable cell shape, cell stacking and incomplete cell overlap while preserving the advantages of LGCA; namely, synchronous transport and binary representation of cells within channels (e.g., a ‘0’ indicating an unoccupied channel and a ‘1’ indicating an occupied channel).

We represent the cells as (1) a single node which corresponds to the position of the cell’s center (or “center of mass”) in the $`xy`$ plane, (2) the choice of occupied channel at the cell’s position designating the cell’s orientation and (3) a local neighborhood defining the physical size and shape of the cell with associated interaction neighborhoods (Figure 3). The interaction neighborhoods depend on the dynamics of the model and need not exactly overlap the cell shape. In our models for rippling and aggregation, we define the size and shape of the cell as a $`3\times \mathrm{}`$ rectangle, where $`\mathrm{}`$ is cell length. As $`\mathrm{}`$ increases, the cell shape becomes more elongated. A cell length of $`\mathrm{}=30`$ corresponds to the $`1\times 10`$ proportions of rippling Myxococcus xanthus cells . Representing a cell as an oriented point with an associated cell shape is computationally efficient, yet approximates continuum dynamics more closely than assuming point-like cells, since elongated cells may overlap in many ways. We have also solved the cell stacking problem, since overlapping cell shapes correspond to cells stacked on top of each other. This cell representation conveniently extends to changing cell dimensions and the more complex interactions of fruiting body formation.

### 3.2 LGCA model for rippling

We assume precise reflection and investigates the roles of a cell refractory period, a minimum response period, a maximum oscillation period and non-linear dependence of reversals on C-factor independently.

#### 3.2.1 Local Rules

1. Our model employs a square lattice with periodic boundary conditions imposed at all four edges. Unit velocities are allowed in the positive and negative $`x`$ directions. (A resting channel may be easily added to model a small percentage of resting cells as in .
2. Cells are initially randomly distributed with density $`\delta `$, where $`\delta `$ is the total cell area divided by total lattice area.
3. Every cell is initially equipped with an internal timer by randomly assigning it a clock value between 1 and a maximum clock value $`\tau `$. We define a refractory period $`R`$ such that $`0R<\tau `$ (see a detailed description of the internal timer, below). If the internal timer $`\varphi `$ of a cell is less than $`R`$, the cell is refractory. Otherwise, the cell is sensitive.
4. At each time-step, the internal timer of each refractory cell is increased by 1 while the internal timer of sensitive cells is increased by an amount proportional to the number of head-on cell-cell collisions $`n`$ occurring at that timestep.
5. When a cell’s internal timer has increased past $`\tau `$, the cell reverses, the internal timer resets to 0, and the cell becomes refractory. Reversals occur as a cell’s center switches from a right- or left-directed channel to a left- or right-directed channel, respectively.
6. During the final transport step, all cells move synchronously one node in the direction of their velocity by updating the positions of their centers. Separate velocity states at each node ensure that more than one cell never occupies a single channel.

#### 3.2.2 Internal timer

We model an internal timer with three parameters; $`R`$, $`t`$ and $`\tau `$. $`R`$ is the number of refractory time-steps, $`t`$ is the minimum number of time-steps until a reversal and $`\tau `$ is the maximum number of time-steps until a reversal. The minimum period of time required for a sensitive cell to become stimulated to turn is the minimum response period $`tR`$. During the refractory period, cells are insensitive to collisions and the internal timer advances at a uniform rate. After the refractory period, cells become sensitive and during this phase the number of head-on cell-cell collisions accelerates the internal timer so that the interval between reversals shortens. This acceleration is density-dependent, so that many simultaneous collisions accelerate the internal timer more than only one collision.

Our internal timer extends the timer in Igoshin et al. . They used a phase variable $`\varphi `$ to model an oscillating cycle of movement in one direction followed by a reversal and movement in the opposite direction. During the refractory period the phase variable advances at a constant rate but during the sensitive period, the phase variable advance may increase non-linearly with the number of collisions. Thus, the evolution of our timer determines reversal rather than a collision as in the model of Börner et al. . The state of our internal timer is specified by $`0\varphi (t)\tau `$. $`\varphi `$ progresses at a fixed rate of one unit per time-step for $`R`$ refractory time-steps, and then progresses at a rate $`\omega `$ that depends non-linearly on the number of collisions $`n`$ which have occurred at that timestep to the power $`p`$:

$$\omega (x,\varphi ,n,q)=1+\left(\frac{\tau t}{tR}\right)\left(\frac{[min(n,q)]^p}{q^p}\right)F(\varphi ),$$
(1)
where,

$$F(\varphi )=\{\begin{array}{cc}0,\hfill & \text{for }0\varphi R\text{;}\hfill \\ 0,\hfill & \text{for }\pi \varphi (\pi +R)\text{;}\hfill \\ 1,\hfill & \text{otherwise.}\hfill \end{array}$$
(2)
This equation is the simplest that produces a reversal period of $`\tau `$ when no collisions occur, a refractory period of $`R`$ time-steps in which the phase velocity is one, and a minimum reversal period of $`t`$ when a threshold (quorum) number $`q`$ of collisions occurs at every sensitive time-step. There is “quorum sensing” in that the clock velocity is maximal whenever the number of collisions at a time-step exceeds the quorum value $`q`$. A particle will oscillate with the minimum reversal period only if it reaches a threshold number of collisions during each non-refractory time-step (for ($`tR`$) time-steps). If the collision rate is below the threshold, the clock phase velocity is less than maximal. However, as the number of collisions increases from 0 to $`q`$, the phase velocity increases non-linearly as $`q`$ to the power $`p`$.

While in the model of Börner et al.  there is no minimum response period for a cell to reverse, and in the model of Igoshin et al.  a minimum response time is an inherent component of the internal clock, our model incorporates “on-off switches” for a refractory period, minimum response period, maximum oscillation period and quorum sensing. Setting the refractory period equal to 0 time-steps in our model is the off-switch for the refractory period, and setting $`t=R+1`$ is the off-switch for the minimum response time. No maximum oscillation period is modeled by choosing a maximum oscillation period $`\tau `$ greater than the running period of the simulation, so that the automatic reversal of cells within $`\tau `$ time-steps has no effect on the dynamics of the simulation. There is no quorum sensing if $`q`$ is set to 1 so that a single collision during a timestep has the same effect as many collisions.

If there is no refractory period, cells are always sensitive to collisions. If there is no minimum response time, cells may reverse immediately after becoming sensitive if there are sufficiently many collisions in one timestep. Finally, if there is no maximum oscillation period, cells may never reverse without sufficiently many collisions.

#### 3.2.3 Head-on cell-cell collisions

We define an interaction neighborhood of eight nodes for the exchange of C-factor at the poles of a cell of length l (see Figure 3). The cell width of 3 nodes is larger than 1 to account for coupling in the y-direction and the interaction neighborhood must extend at least two nodes along the length of the cell to compensate for the discretization of the lattice since cells traveling in opposite directions may pass without their poles exactly overlapping.

A head-on cell-cell collision is defined to occur when the interaction neighborhoods of two anti-parallel cells overlap. A cell may collide with multiple cells simultaneously since the interaction neighborhood is four nodes at each pole. Note that the specific shape of the cell is not important for rippling dynamics since the two areas of C-signaling are the only places where interaction occurs. Nevertheless, a shape extending over several nodes is necessary to permit the necessary overlapping and stacking at high density since the exclusion principle mandates that each channel has at most one cell center. Thus, the cell centers of two colliding cells will be separated by one cell length and do not compete for channels at the same node. Also, for sufficiently long cell lengths, the probability of more than one cell center located at the same node is low even when the local cell density may be high.

We are able to simulate a rippling population with arbitrary concentrations of both wild-type and non-C-signaling cells and quantitatively reproduce their experimental results in detail, as did Igoshin et al. using their continuum model . Further, we demonstrate that the change in wavelength may be entirely explained by the change in density of C-signaling cells.

## 4 Rippling results and discussion

Our model forms a stable ripple pattern from a homogeneous initial distribution for a wide range of parameters, with the ripples apparently differing only in ripple wavelength, ripple density and ripple width (see Figure 4).

Absence of a maximum oscillation period is modeled by choosing a maximum oscillation period $`\tau `$ greater than the running period of the simulation, so that the automatic reversal of cells within $`\tau `$ time-steps has no effect on the dynamics of the simulation. We find that ripples form with or without a maximum oscillation period over the full range of densities. When there is a maximum oscillation period, the maximum oscillation period must be chosen greater than twice the refractory period for the development of ripples. There is no upper bound on the maximum oscillation time, which is why the maximum reversal period is unnecessary. Ripples develop most quickly and cell oscillations are most regular with an internal timer when the maximum oscillation period is carefully chosen with respect to the other parameters of the model. Nevertheless, it appears that experimental results are best reproduced when there is no maximum oscillation period.

A refractory period is required for rippling for cells of length greater than 2 or 3 nodes, and although there may exist a minimum response time of more than one time-step, it is an interesting result of our model that the minimum response period $`TR`$ must be small compared to the refractory period. In particular, rippling occurs whenever the minimum oscillation time $`t`$ is greater than $`\mathrm{}/v`$ time-steps and the refractory period $`R`$ is at least two-thirds $`t`$. The first condition is required because if the minimum oscillation period $`t`$ is less than the period of time it takes a cell to travel one cell length, two cells or a cluster of cells will stimulate each other to oscillate in place. The second condition that the refractory period is at least two-thirds the minimum oscillation period indicates that the minimum response time of a cell can not be too long compared to the refractory period.

Experiments suggest that the minimum oscillation period of a cell in response to C-factor is about 3 minutes . According to our result that the minimum response time cannot be more than two-thirds the refractory period, we can predict the existence of a refractory period in myxobacteria cells, with a duration of 2-3 minutes.

The wavelength of the ripples depends on both the duration of the refractory period and the density of signaling cells. Figure 5 shows that the ripple wavelength increases with increasing refractory period (a) and decreases with increasing cell density (b). Notice that error bars that show standard deviations of the mean wavelength over five simulations increases with wavelength. A refractory period of 2-3 minutes yields a ripple wavelength of about 60 micrometers (Figure 5a), which corresponds well to typical experimental ripple wavelengths . The correspondence between refractory period and wavelength given in Figure 5 is a only rough estimate, however. We believe the reasons are that in these simulations the cell density is relatively low, which decreases the density of C-factor relative to experimental conditions, and cells are not very elongated, which increases the density of C-factor relative to experimental conditions.

Note that in Figure 5a, the curve has a wavelength of approximately 20 micrometers when the refractory period is less than 1 minute. Since cells have a length of 5 micrometers, this is the smallest wavelength that may be resolved as there is only one cell length between subsequent high density waves. At very high density, when the refractory period is 0, cells may be stimulated to reverse every timestep, so that there would be, theoretically, a wavelength of only 1 node. However, cells will be uniformly distributed in this case and there will be no well-defined high-density waves. In the simulations described in Figure 5b, density is increased while refractory and minimum oscillation periods have a constant value. The minimum possible wavelength in this case is limited by the minimum oscillation period. In particular, the minimum possible wavelength is twice the minimum distance traveled by a cell between reversals, which is twice the distance traveled during the minimum oscillation period, which is 30 micrometers in this example. Thus, even as density is increased very high, the curve must have a horizontal asymptote at wavelength = 30 micrometers.

### 4.1 Non-linear response of reversals to C-factor.

Reversals depend on the number of collisions a cell encounters which depends on the density of C-factor. Thus the number of collisions required for a reversal, the quorum value $`q`$, should be a function of the density of C-signaling nodes. The density of C-signaling nodes is a function of both cell density and cell length since longer cells have a reduced C-signaling area to non-C-signaling area ratio. Thus, we describe optimal quorum values $`q`$ as a function of C-signaling node density rather than cell density.

At a low density of C-signaling nodes, ripples form even when both $`q`$ and $`p`$ are 1 so that only 1 collision during the sensitive period is needed to trigger an reversal. When the density of C-signaling nodes is greater than or equal to 1, however, the chances of collisions are so high in the initial homogeneous population that cells almost always reverse in the minimum number of timesteps and, with no differential behavior among cells, a rippling pattern fails to form. Ripples will form at arbitrarily large densities of cells and C-signaling nodes if the number of collisions needed to trigger a reversal is increased. When the number of collisions required for a reversal is greater than 3 ($`q>3`$), ripples develop more quickly if the non-linear response to density $`p`$ is increased greater than 1. A value of $`p=3`$ yields optimal rippling for all quorum values and densities, which is consistent with the results of Igoshin et al.  for their value of $`p`$ .

### 4.2 Ripple phase shift

Counter-propagating ripples appear to pass through each other with no interference, which lead Sager and Kaiser to propose the hypothesis of precise reflection . Indeed, tracking of right-propagating ripples and left-propagating ripples in Figure 8a, shows that the waves move continuously despite collisions and subsequent reflection. Inspection of the collision and subsequent reversal of two cells, however, shows there is a jump in phase equal to exactly one cell length if they reverse immediately upon colliding (see Figure 6a). This phase jump occurs because a cell reverses by changing its orientation rather than by turning: when a right-moving particle collides with a left-moving particle and reverses, it is exactly one cell length ahead of the left-moving cell that it replaces. When all of the particles within a ripple are in phase, as is often the case, this jump is also seen in the ripple waves as two waves interpenetrate. If the cells continue $`p`$ more steps before reversing (for example, if their clocks were almost near $`\tau `$ after the collision), then there would be a phase jump of $`\mathrm{}2p`$. If $`2p>\mathrm{}`$, there will be a phase delay (see Figure 6b). In their continuous model, Igoshin et al. (, Figure 3b) also showed when ripples collide a small jump in phase reminiscent of a soliton jump.

### 4.3 Effect of dilution with non-signaling cells

Sager and Kaiser  diluted C-signaling (wild-type) cells with non-signaling (csgA minus) cells that were able to respond to C-factor but not produce it themselves. When a collision occurs between a signaling and a non-signaling cell, the non-signaling cell perceives C-factor (and the collision), whereas the C-signaling cell does not receive C-factor and behaves as though it has not collided. The ripple wavelength increases with increasing dilution by non-C-signaling cells. Simulations of this experiment with and without an internal timer with a maximum oscillation period give very different results. Figure 7a shows that the dependence of wavelength on the fraction of wild type cells resembles the experimental curve only when there is no maximum oscillation period assumed in the model (compare with Figure 7.G in .) Thus our model predicts that rippling cells do not ripple with a maximum oscillation period. Notice that the range of wavelengths when there is no assumed maximum oscillation period is in good quantitative agreement with that of experiment (compare Figure 7a, solid line with Figure 7.G in .

Igoshin et al.  have previously reproduced the experimental relationship between wavelength and dilution with non-signaling cells (see  Supplemental materials, Figure 8) by adjusting their original internal timer. As the density of C-signal decreases, the phase velocity slows linearly and the maximum oscillation period of the internal timer increases continuously. Thus, the maximum oscillation period varies in their model. We assume a constant maximum oscillation period, which is either present or absent (longer than the simulation running time). Note that if the maximum oscillation period increases sufficiently with decreased density of C-factor so that a cell is always stimulated to turn before the internal timer would regulate a turn, then the addition of an internal timer is superfluous. In this case, the two models are similar.

Our simulations show ripple wavelength increases with increased dilution by non-signaling cells. Since wavelength also increases with decreasing density of signaling cells (Figure 5b), we ask if the mutant cells have any effect on the rippling pattern. Figure 7b shows the wavelength dependence on the density of signaling cells when only signaling cells are present (dotted line) and for a mixed population of signaling cells of the same density with non-signaling cells added so that the total cell density is always $`1.6`$ (solid line). Apparently, the decrease in C-factor explains the increase in wavelength. The non-signaling mutants do not affect the pattern at all.

As density increases, wavelength decreases and the larger number of cells are distributed over a greater number of ripples. This result is further evidence of the the role C-signaling plays as a density-sensing and density-modulating mechanism. To test this further, we ran a simulation for initial conditions in which a high density stripe stretches vertically down a lattice. As ripples formed and propagated, the cells were quickly distributed more evenly over the lattice (Figure 8b). The redistribution of cells occurs much faster than if cells moved randomly at each time-step (compare Figure 8, b and c). Thus, although there is no net transport of cells larger than one wavelength when cells are evenly distributed , there is net migration of cells away from high density regions.

### 4.4 Discussion

In our simulations, high-density waves of cells form from a homogeneous distribution of cells for a wide range of parameters and initial cell density. At low cell density when there is no assumed maximum oscillation period, the wavelength between ripples is large. The explanation for this is that a larger region of the lattice must be “swept” to collect an aggregate of cells with sufficiently high density to reverse another aggregate. In the extreme case where density is nearly zero, a single cell will keep traveling without ever encountering sufficient collisions, and the wavelength is infinite. Thus, the mechanism of rippling may be viewed as a mechanistic “sweeping” of an arbitrarily large area, in which cells modulate the area that they span between reversals so as to efficiently collect ripples of a minimum density.

At very high density, by the same argument, wavelength is small. The number of ripples per area increases so that the number of cells per ripple does not increase linearly with the increase in density, but less. Nevertheless, ripples formed from high density initial conditions do result in ripples which are wider and more dense than ripples formed from low density initial conditions. This may not be viewed as a limitation in design, but as the foundation of another possible role for rippling in myxobacteria: the even distribution of cells over an arbitrarily large area. When there are both high and low cell density regions (as in the simulation of Figure 8b), many high-density waves will form in the region of highest density, and fewer, lower density waves will form in the regions of lower density. At the interface between these regions, a high density wave encounters a lower density wave. The high density wave will reverse all the cells of the low density wave such that all the cells from the low density wave return toward the low density region. The cells of the low density wave, on the other hand, will not be able to reverse all the cells of the higher density wave. Rather, the lower density wave will only be able to reverse a proportional number of cells in the high density wave. The surplus cells in the high density wave will continue without reversing into the low density region. Thus, the rippling mechanism in a region of variable cell density creates a “pulse” of surplus cells within the high density region which is efficiently directed into lower density regions.

In summary, rippling is an efficient mechanism for both forming evenly spaced accumulations of high cell density, and evenly spaced accumulations of nearly equal cell density. In experiments, cells do not reflect by exactly 180 degrees. However, since most cells move roughly parallel to each other, models based on reflection are reasonable approximations. Modeling the experimental range of cell orientations would require a more sophisticated CA since LGCA require a regular lattice which does not permit many angles. In the aggregation section below we describe a model on a triangular lattice which could be adjusted to incorporate 120 degree reversals. Although a cell is ready to turn when the internal timer $`\varphi `$ is greater than $`\tau `$, it may not be able to turn if the opposite channel is already occupied. This is another limitation of our LGCA model. We handle this situation by continuing to transport the cell in its direction of orientation at each time-step until the opposite channel is available. The effect of this delay is negligible, even at high densities within a ripple, when cells are so long that the probability of two cell centers occupying the same node is very small.

## 5 A preliminary model for aggregation

Rippling is an intermediate, transient stage of fruiting body formation, which is not necessary for aggregation formation . Figure 1 shows a field of aggregation centers surrounded by ripples. In this section we present a different LGCA model based on C-signaling alignment. This LGCA model reproduces the sequence and geometry of the non-rippling stages of fruiting body formation in detail, demonstrating how C-signaling-based alignment can account for these patterns with very few additional assumptions.

The non-rippling stages of fruiting body formation include alignment, streaming and aggregation. During alignment, cells form aligned patches from a random distribution. While streaming, cells form long chains which move cooperatively into aggregation center . Aggregation is the phase in which cells form rounded collections that may either recede or mature into fruiting bodies. We model aggregation including only a simple local rule for C-signal-based alignment. The aggregates formed in our model are not species-specific and do not include local rules for rippling.

We use a hexagonal lattice since cell motion during aggregation is not one-dimensional as in rippling. In this specialized LGCA model, identical rod-shaped cells are all modeled as $`3\times \mathrm{}`$ rectangles with C-signal interaction neighborhoods as depicted in Figure 9. Cells move exactly one node per time-step in the direction of their orientation and there may only be one cell center per channel per node. In contrast to rippling, we find that the cell aspect ratio is an important parameter for streaming and aggregation, the simulations presented here all have a $`1\times 7`$ aspect ratio for cells.

Myxobacteria align when they move. We choose an alignment based on C-signaling. We use a Monte Carlo process, in which cells turn 60 degrees clockwise, 60 degrees counter-clockwise or persist in their original direction with probability favoring directions that maximize overlap of C-signaling nodes. Interaction only occurs when the C-signaling nodes at the head of a cell overlap with the C-signaling nodes at the tail of another cell, and interaction occurs regardless of cell orientation. Head and tail C-signaling neighborhoods are shown in Figure 9.

We model C-signal alignment on a $`256\times 256`$ lattice, in which our initial conditions are a random distribution of cells at high density. Within a few time-steps, cells form aligned patches (see Figure 10a). This reproduces the initial alignment stage of myxobacteria cells during fruiting body formation in which cells form an aligned patchwork . Since there are only 6 directions permitted on the lattice, the aligned patchwork appears as a very regular triangular network. Cells are aligned both parallel and anti-parallel within each patch since the overlap of C-signaling nodes is maximized when cells are aligned regardless of their orientation. This geometry is significant because cells have naturally formed the aligned bi-directional arrays necessary for rippling. The local rules for rippling have been suppressed in this preliminary model, however, to evaluate the patterns formed by C-signaling-based alignment alone.

The homogeneous triangular network is not stable over time. As cells move, they turn and flow from one patch to another. Cells moving in a low density area are likely to turn into a higher density patch, so the network of cells condense into thick aggregates (Figure 10b). In experiments, a myxobacteria stream will often merge into an adjacent myxobacteria stream . Figure 10b shows that the cells move into the aggregates along streams directed into the aggregate, reproducing the streaming stage in which myxobacteria cells form long aggregates that move cooperatively .

The aggregates continue to condense while arranging and dividing into many small, circular orbits about 1.5 cell lengths (about 10-20 $`\mu `$m) in diameter (Figure 10c). The aggregates often form in clusters of two or three closed orbits (Figure 10b), corresponding to fused aggregates (sporangioles). In Stigmatella erecta, several fruiting bodies may form in groups and fuse . A magnified picture of the cell centers of a typical aggregate show that cells are arranged in dense, concentric layers tangent to a relatively low-density inner region (Figure 10d). Thus, they are geometrically equivalent to the basal region of aggregates in Myxococcus xanthus.

Cells in our simulation simultaneously move both clock-wise and counter-clockwise around the aggregate, as they do in the fruiting bodies of Myxococcus xanthus . The microscopic circular or elliptic orbits of Stigmatella spp. often disappear as cells spiral away from the aggregate . Similarly most orbits in our simulation also eventually disappear as cells spiral away from the aggregate. Nevertheless, orbits typically survive for several hundred time-steps, which is about 5 to 10 complete revolutions for each cell.

During myxobacteria aggregate formation, several aggregation centers will form and, inexplicably, one aggregation center will grow as an adjacent aggregation center disappears . Our simulations offer a closer look at this process: a stream may form that connects two adjacent aggregation centers and, stochastically, cells stream from one aggregation center to another until the largest aggregation center absorbs the smaller one (Figure 10e).

Figure 10f shows several stable aggregates which have developed at 200 time-steps. Notice that the stable hollow aggregate has a much thicker annulus of cells than the non-stable orbits of Figure 10c, suggesting that only large hollow aggregates are stable. In our simulation at a threshold density within the aggregate, the hollow region of an aggregate will fill with cells such that cells are are arranged in six dense overlapping layers (Figure 10f). The third aggregate shown is a chain of cells which span the lattice and thus form an orbit due to the periodic boundary conditions of the lattice.

It has been proposed in  that circular motility at aggregation sites, trail following and local accumulation of slime account for fruiting body formation in Stigmatella spp. We hypothesize that once C-signaling has drawn cells into an aggregate, cell and slime cohesivity cause myxobacteria cells to round up into a mound while constant cell velocity pushes cells toward the surface of the mound, so that cells form a dense hemisphere of cells spiraling around a hollow center. In our model, a closed, circulating orbit of cells is the only stable configuration of a stationary aggregate since cells are constantly moving. At low and intermediate density, these orbits are hollow and cells are arranged tangentially within an annulus. At a threshold density, however, every channel of every node within the aggregate becomes occupied and there is no hollow center. We hypothesize that in a more advanced three-dimensional model, the addition of a local rule accounting cell and slime cohesivity will cause cells to round while maintaining the hollow center.

This model is preliminary because only C-signal based alignment is modeled and the aggregates formed are not species-specific. In this simulation, a 256x256 lattice size was chosen, which corresponds to an 80x80 $`\mu `$m region and about 10,000 cells for an averaged cell density of 10, much smaller compared to normal fruiting bodies that may be up to 1000 $`\mu `$m in diameter. Aggregation in this model is described in more detail in .

This model can be adjusted to model rippling and aggregation concurrently by incorporating the local rules for rippling described in Section 3.2. Tracking of rippling cells in experiment (see Figure 6 in  and Figure 6 in ) suggest that reversals are about 150 degree changes in orientation rather than exactly 180 degrees, as we have assumed in the model for rippling of this paper. On a triangular lattice, a reversal of 120 degrees would be an equally good approximation of the 150 degree rotation. C-signaling based alignment combined with local rules for rippling would ensure that the majority of cells still remain parallel. Nevertheless, a regular shift in orientation by 120 or 180 degrees may have an interesting effect on the final distribution of cells and, subsequently, fruiting bodies.

Fruiting bodies among different myxobacteria species are very diverse (see Figure 3 in ). For example, while Myxococcus fruiting bodies are a relatively simple, single mound of cells, other species form clusters of mounds called sporangioles that are raised on a stalk. The interaction of adventurous verses social motility may account for these different morphologies . During stalk formation of Stigmatella spp., cells are arranged perpendicularly to the mound as the fruiting body is lifted . Also, fruiting body stalks may be composed of a larger, second cell-type . Thus, once the stages of fruiting body formation have been modeled in general, it would be interesting to determine which parameters need to be varied to model the fruiting bodies of different species.

## 6 Summary

In this paper, we present a new LGCA approach for modeling cells which is computationally efficient yet approximates continuum dynamics more closely than assuming point-like cells. As an example of this new approach, we present a model for myxobacteria rippling based on the hypothesis of precise reflection. The results of our model show that rippling is stable for a wide range of parameters, C-signaling plays an important role in modulating cell density during rippling, and non-C-signaling cells have no effect on the rippling pattern when mixed with wild-type cells. Further, by comparing model results with that of experiment, we can conclude reversals during rippling would not be regulated by a built-in maximum oscillation period. We also present a second LGCA model based on C-signal alignment which reproduces the sequence and geometry of the non-rippling stages of fruiting body formation in detail, showing that a simple local rule based on C-signaling can account for many experimental observations.

Acknowledgments

We would like to thank Dale Kaiser, Frithjof Lutscher and Stan Marée for very helpful discussions. MSA is partially supported by grant NSF IBN-0083653. YJ is supported by DOE under contract W-7405-ENG-36. MAK is supported by the Center for Applied Mathematics and the Interdisciplinary Center for the Study of Biocomplexity, University of Notre Dame, and by DOE under contract W-7405-ENG-36.