text:
# A Search for Short Time-Scale Optical Variability in the GRB 030329 Afterglow

## 1 Introduction

The recent detection of early fluctuations about the mean power-law decay of GRB afterglows Bersier et al. (2002); Mirabal et al. (2003) may provide important clues in our understanding of the afterglow evolution. However, since several interpretations have been put forth to explain the fluctuations (e.g. Heyl & Perna (2003)), it is crucial to determine the mechanism ultimately responsible for the variability. Here we present an initial attempt to help constrain the microphysics of GRB afterglows using a rapid optical photometric sequence.

## 2 Discussion

High-speed optical photometry of the GRB 030329 afterglow was obtained at the MDM 1.3 m telescope beginning on March 30 03:05 UT, 15.5 hours after the burst. The data consists of a 7.5–hour time series in the $`R`$-band with a time resolution of 90 s. Altogether, we obtained a total of 306 points. During this period, the magnitude of the optical afterglow declined from $`R=15.4`$ to $`R=16.2`$. Figure 1 shows a power-law decay slope $`\alpha `$ = –1.931 $`\pm `$ 0.005 fitted to the data, where the uncertainties are dominated by faint comparison stars and large airmasses at the beginning and the end of the run.

The computed residuals indicate that the average fluctuation about the power-law fit is smaller than 0.01 mag on time scales of minutes to hours. If the flux variability $`f_1/f_0`$ traces the density contrast $`n_1/n_0`$ as a function of $`f_1/f_0(n_1/n_0)^{1/2}`$ Lazzati et al. (2002), this implies that the enhancements in the vicinity of the GRB are less than a factor of 1.03 in density. In other words, the circumburst medium was close to homogeneous between 1.4 $`\times 10^{17}`$ cm and 1.7 $`\times 10^{17}`$ cm from the burst site, assuming a spherical adiabatic blast wave expanding in a stellar wind. Unfortunately this fast sequence by itself cannot place strong limits on intrinsic afterglow variability, but it shows that the afterglow evolution can be highly uniform on short time-scales.

## 3 Conclusions and Future Work

We have shown preliminary results of a high-speed optical photometry program to access short time-scale fluctuations in GRB afterglows. A comprehensive effort to observe and model fluctuations at different time scales might hold the key to understanding the microphysics of GRB afterglows.


format: text

text:
# Spontaneous violation of CP symmetry in the strong interactions

## Acknowledgements

This manuscript has been authored under contract number DE-AC02-98CH10886 with the U.S. Department of Energy. Accordingly, the U.S. Government retains a non-exclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for U.S. Government purposes.


text:
# Coupled Folding-Binding versus Docking: A Lattice Model Study

## I Introduction

Protein structures are often viewed as well-defined static entities. For many proteins, this simplified static picture is accurate enough to provide valuable information about the function.Branden:91  However, there are proteins that are wholly unstructured or only partially structured and yet functional. In fact, recent studies suggest that such proteins are more common than previously thought, especially in eukaryotic cells;Wright:99 ; Dyson:02  for instance, it has been estimated that as much as 17% of the proteins in Drosophila are wholly unstructured.Dunker:01  In many cases, intrinsically unstructured proteins adopt specific structures upon binding to their biological targets.Wright:99 ; Dyson:02  Folding and binding are then coupled, thus establishing a direct link between folding and function.

It has been suggested that unstructured proteins offer several advantages in cellular regulation.Wright:99 ; Dyson:02  Being unstructured might, for instance, allow one protein to interact with several targets. It might also be useful for control purposes due to rapid turnover. Furthermore, it has been argued that being unstructured might facilitate the binding of the protein to a target,Shoemaker:00  by increasing the ‘capture radius’. This mechanism, termed ‘fly casting’,Shoemaker:00  was analyzed using a low-dimensional representation of the binding process. More recently, the thermodynamics of coupled folding-binding were examined using a generalized random energy model.Papoian:03

Here we study coupled folding-binding by computer simulations of a simple chain-based model. Simulating coupled folding-binding of a chain in a controlled manner requires proper sampling of the full conformational space and is harder than simulating folding of an isolated chain, since rigid-body translations and rotations must be taken into account. The complexity of the problem makes it highly desirable to study coarse-grained models before entering high-resolution modeling. For the present study, we use the minimal two-dimensional hydrophobic/polar HP lattice model of Lau and Dill.Lau:89  This model has been widely used to investigate basics of protein folding.Chan:02  This and similar lattice-based models have also been used to gain insights into topics such as protein evolution,Bornberg-Bauer:97 ; Cui:02 ; Blackburne:03  prion-like conformational propagationHarrison:99 ; Harrison:01  and protein aggregation.Broglia:98 ; Istrail:99

Our study is inspired by recent experiments by Wahlberg et al.Wahlberg:03  on the in vitro evolved protein ZSPA-1. This sequence was engineeredEklund:02  from the Z domain of staphylococcal protein A, a well characterized three-helix-bundle protein.Tashiro:97  It was selected for binding to the Z domain itself. In the Z:ZSPA-1 complex, ZSPA-1 adopts a structure similar to the solution structure of the Z domain.Wahlberg:03 ; Hogbom:03  However, in solution, ZSPA-1 turns out to be structurally disordered.Wahlberg:03  The engineered ZSPA-1 thus exhibits coupled folding-binding. In Ref. Favrin:03, , computer simulations of the solution behaviors of ZSPA-1 and the Z domain were performed, using a relatively detailed off-lattice model with 5–6 atoms per amino acid.Favrin:03  Simulating the binding behavior of an unstructured protein like ZSPA-1 at this level of resolution remains, however, a challenge.

The present study consists of two parts. First, we study binding statistics for very short HP chains with up to $`N=14`$ monomers, to get an idea of how likely binding is to occur for stable and unstable sequences, respectively, in this model. We then perform a more detailed study of two $`N=25`$ sequences. In particular, this study allows us to compare the free energy of coupled folding-binding with free-energy profiles for docking and for folding of an isolated chain.

## II Model and Methods

In the HP model,Lau:89  the protein chain is represented by a string of hydrophobic (H) or polar (P) beads on the square lattice. Adjacent beads along the chain are connected by links of unit length. It is forbidden for two beads to occupy the same lattice site. Two beads that are neighbors on the lattice but non-adjacent along the chain are said to be in contact. The energy that a configuration gets is determined by the number of HH contacts, each HH contact being assigned an energy $`ϵ<0`$. This defines the model for a single chain in isolation.

Here we study single chains interacting with some fixed target structure. Two types of targets are considered. In Sec. III.1, the target is an immobilized HP chain, and the inter-chain interactions are taken to be the same as the intra-chain interactions; that is, each HH contact between the target and the flexible chain is assigned the energy $`ϵ`$. Any cross HH contact is given this energy whether or not it is present in the final bound state. In Sec. III.2, the target is an extended structure with one particular binding site. In this case, there is only one specific bead of the flexible chain that can interact with the binding site of the target. A contact between these two sites is taken to be favorable by an energy $`ϵ_\text{b}=3ϵ`$.

To study the thermodynamics of these systems, we use Monte Carlo methods. The moving chain is allowed four types of moves: local one- and two-bead moves, non-local pivot moves, and one-step translation moves (left, right, up or down). All moves are subject to a Metropolis accept/reject step, to ensure that detailed balance is fulfilled. One round of simulation, called one ‘sweep’, consists of (at most) $`N1`$ one-bead moves, $`N2`$ two-bead moves, one pivot move and one translation move, $`N`$ being the number of beads.

In Sec. III.2, in addition to the thermodynamic simulations, we also study the binding process as a function of Monte Carlo time. These calculations follow exactly the same protocol, except that the non-local pivot update is omitted, to avoid large unphysical deformations of the chain.

With its simplified conformational space and its minimal two-letter alphabet, the HP model is not meant for studies of specific proteins, but rather to shed light on general questions for generic sequences. It is worth notingChan:02  that if the principle of minimal frustrationBryngelson:87  holds, then the native structures of functional proteins should be strongly favored by the hydrophobicity pattern alone. This suggests that the HP model, despite its simplicity, might be able to capture non-trivial features of the mapping from sequence to native structure. A statistical analysis of HP model sequences with unique ground states lends support to this view; it turns out that the hydrophobic beads are anticorrelated along the chains,Irback:00 ; Irback:02  which is the same behavior that real (globular) protein sequences show.Irback:96 ; Irback:00  A recent study of the distribution of hydrophobicity in protein sequences can be found in Ref. Sandelin:03, . For a study of designed hydrophobic/polar copolymers with positive hydrophobicity correlations, see Refs. Khokhlov:98, ; Govorun:01, .

The HP model has the advantage over high-resolution models that exact results can be obtained for short chains, by exhaustive enumeration of all possible states. On the square lattice, all possible sequences with unique ground states (in isolation) have been determined for $`N25`$,Irback:02  along with the corresponding structures. Below, we make use of these results. Sequences having a unique ground state will be referred to as stable.

## III Results and Discussion

### III.1 Statistics for Short Chains

To get an idea of how likely different binding behaviors are to occur in this model, we start by studying a large number of short sequences with up to $`N=14`$ beads. Each sequence interacts with some immobilized target chain with the same length. All target sequences considered are stable and held fixed in their unique ground-state conformations. The two-chain system is contained in a box of size $`2N\times 2N`$. The target is fixed in the middle of the box so that the moving chain can attack it from any side. We say that the moving protein binds to the target if there is a unique minimum-energy configuration for the two-chain system. Whether or not this criterion is met can be determined with high confidence by Monte Carlo methods for these chain lengths. As an operational criterion for binding we require that the minimum-energy configuration is unique and is visited at least 10 times during the course of the simulation. To ensure that these visits to the bound state are ‘independent’, a visit is counted only when the system comes to this state after going to some state which had at least 20% higher energy than the bound state.

We first test all possible sequences with a given $`N`$ for binding to one particular target conformation with that $`N`$, using $`N=10`$, 11 and 12. As target, we somewhat arbitrarily pick the first entry in the list of all stable ground states obtained in Ref. Irback:02, .footnote  The results of these calculations are summarized in Table 1. From this table it can be seen that the fraction of binding sequences is a few per cent, which is comparable to the number of stable sequences. A clear majority of the binding sequences are unstable, which means that docking cannot be the mode of binding. That most binding sequences are unstable does not mean that such sequences have an intrinsically higher propensity to bind than stable ones, but rather it merely reflects the fact that most sequences are unstable.

As mentioned in the introduction, one suggested use of being unstructured is that it might enable the chain to bind to different targets. In the model, there are sequences showing this behavior. An example of this can be found in Fig. 1, which shows an unstructured $`N=12`$ sequence that is able to bind to two different targets (both of which are held fixed). The structure that the unstable sequence adopts upon binding is seen to depend on the target.

Finally, restricting ourselves to stable sequences, we test for self-binding, by using the same sequence as target and moving chain simultaneously. In obtaining Table 1 above, this calculation was done for those three sequences that served as targets. Using $`N=12`$ and $`N=14`$, we now consider all stable sequences. Table 2 shows the results of these calculations. We see that $``$10% of the stable sequences are self-binding both for $`N=12`$ and $`N=14`$. For $`N=12`$, we observe that all the 9 self-binding stable sequences bind in their ground-state conformations. For $`N=14`$, on the other hand, there exist 2 self-binding stable sequences for which the bound structure differs from the isolated one, while the remaining 40 keep the same structure after binding.

Let us briefly summarize the results presented so far. Our study of all possible sequences for $`N=10`$, 11 and 12 shows that a significant fraction (a few per cent) can bind to the targets considered. Furthermore, most of these binding sequences are unstructured in isolation. Our study of all stable sequences for $`N=12`$ and $`N=14`$ shows that $``$10% of them are self-binding. For most of these self-binding sequences, the complexed conformation is identical to the stable structure in isolation, indicating the possibility of a docking-like binding behavior.

### III.2 Coupled Folding-Binding versus Docking

Having studied binding propensities of short chains, we now turn to a more detailed study of two longer chains with $`N=25`$, whose binding behaviors are contrasted. This time the target is an extended structure located in the left-bottom end of an $`N\times N`$ box, in which the chain is confined (see Figs. 2 and 3). The target has only one binding site, and the binding is specific in that only one bead of the moving chain can interact with this binding site.

The two sequences considered are given in Table 3. The first sequence, called S, was studied in isolation in Ref. Irback:02, . It was obtained by applying a sequence optimization algorithmIrback:98  to 326 stable sequences all having the same unique ground-state conformation. It turns out that S has a ground-state energy of 13$`ϵ`$. The next most favorable conformations have only 11 HH contacts, so there is an energy gap of $`2ϵ`$, making S unusually stable. The other sequence studied, called U, is unstable. Its degenerate ground state has an energy of $`11ϵ`$. In our calculations, we observed five different conformations of U with this energy. Sequence U is a close analog of sequence S with just one hydrophobic bead mutated to a polar one. The two sequences were taken to be similar in order for their bound structure to be the same.

Our thermodynamic simulations of these two sequences were started from random configurations and contained $`10^8`$ Monte Carlo sweeps each. The results show first of all that both these sequences do bind to the chosen target, in the sense that the system of target plus moving chain has a unique minimum-energy state. This state, referred to as the bound state, was visited many independent times in the simulations. The bound structure is the same for both chains (see Figs. 2 and 3), and coincides with the stable structure of S in isolation. The energy of the bound state is $`13ϵ+ϵ_\text{b}`$ for S and $`11ϵ+ϵ_\text{b}`$ for U.

To characterize the binding behaviors of these sequences, we monitor the following two quantities:

1. The binding parameter $`I`$. $`I`$ is the geometric distance between the binding site of the target and that bead of the moving chain that can interact with the target. As binding progresses, the value of $`I`$ reduces.
2. The folding parameter $`Q^{}`$. The chains studied form 14 internal contacts, called native contacts, in their bound states. The number of these native contacts being present, $`Q`$, provides a measure of the ‘nativeness’ of the chain. As the chain folds to its bound-state structure, the value of $`Q`$ increases. It turns out that the value $`Q=13`$ is impossible to attain for these sequences. In our free-energy calculations, we therefore use a folding parameter $`Q^{}`$ defined by

$$Q^{}=\{\begin{array}{cc}Q\hfill & \mathrm{if}Q12\hfill \\ 13\hfill & \mathrm{if}Q=14\hfill \end{array}$$
(1)

Figures 4 and 5 show the free energy calculated as a function of these two variables, $`F(I,Q^{})`$, for the sequences S and U, respectively. The free energy is defined by $`P(I,Q^{})=\mathrm{exp}[F(I,Q^{})/kT]`$, where $`P(I,Q^{})`$ is the joint $`I,Q^{}`$ probability distribution. The temperature is taken as $`T=ϵ/2.6k`$ for S and $`T=ϵ/3.2k`$ for U ($`k`$ is Boltzmann’s constant). These temperatures are chosen such that the binding probabilities for the two chains are significant and close to each other ($`35`$ %). Note that the systems are in their minimum-energy states if and only if $`(I,Q^{})=(0,13)`$.

For sequence S (Fig. 4),

we find that $`F(I,Q^{})`$ has a simple shape with two narrow valleys along the lines $`Q^{}=13`$ and $`I=0`$, respectively, and a broad and shallow minimum centered at $`(I,Q^{})(11,5)`$, where the chain is unbound and unfolded. This suggests that there are two very different major binding modes for this sequence. One way for the chain to reach its bound state is along the line $`Q^{}=13`$, which corresponds to rigid docking; the chain first folds and then moves towards the binding site. The other major binding mode is along the $`I=0`$ valley. Here the chain first attaches to the binding site and then folds to its final shape. Following Ref. Shoemaker:00, , we refer to this behavior as the fly-casting mechanism.

Let $`P_\text{d}`$ and $`P_\text{f}`$ denote the probabilities of finding the system in the $`Q^{}=13`$ and $`I=0`$ corridors, respectively (not counting the bound state); that is,

$$P_\text{d}=\underset{I>0}{}P(I,Q^{}=13)P_\text{f}=\underset{Q^{}=0}{\overset{12}{}}P(I=0,Q^{})$$
(2)
where the subscripts d and f refer to docking and fly casting. For sequence S, we find that $`P_\text{d}`$ is about twice as large as $`P_\text{f}`$ (see Table 4), suggesting that docking is the preferred binding mode for this sequence.

For the unstructured sequence U (Fig. 5),

there is no free-energy valley corresponding to rigid docking, so $`P_\text{d}`$ is small (see Table 4). There is, by contrast, an $`I=0`$ valley corresponding to fly casting for this sequence, too. The population $`P_\text{f}`$ is, in fact, higher for U than for S (see Table 4). The $`I=0`$ valley is, for both sequences, separated from rest of the conformational space by a free-energy ridge. For sequence U, there are some narrow troughs in these hills towards the $`I=0`$ corridor, which should make it easier for this sequence to reach this corridor. In addition to the $`I=0`$ valley, sequence U exhibits two somewhat less pronounced valleys along the lines $`Q^{}=9`$ and $`Q^{}=12`$, respectively, as well as a broad unbound and unfolded minimum at lower $`Q^{}`$. Note that there exist $`Q^{}=12`$ conformations with minimal intra-chain energy for this sequence. The presence of the $`Q^{}=9`$ and $`Q^{}=12`$ valleys suggests that this chain can follow several different paths to its bound state. If the chain follows one of the fixed-$`Q^{}`$ valleys, it reaches the binding site without having its full bound-state structure; folding is completed after the chain has attached itself to the target. When increasing the box size to $`3N\times 3N`$, we observed valleys at the same values of $`Q^{}`$ which were stretched along $`I`$.

It is instructive to take a closer look at the free energy in the docking and fly-casting corridors, respectively. Figure 6a shows the free energy along the docking corridor, $`F(I,Q^{}=13)`$, for sequence S. We see that in order to reach the bound state, the chain has to pass a free-energy barrier with a height of $``$ 3 $`kT`$. This barrier is entropic. Figure 6b shows the free energy along the fly-casting corridor, $`F(I=0,Q^{})`$, for both sequences. Two features are worth noting. First, the shape of the curve is roughly the same for the two sequences, although the curve is shifted downwards for sequence U which spends more time in this corridor. Second, there is no major free-energy barrier in this corridor (the highest barrier is about 1 $`kT`$), and the last part towards the bound state (from $`Q^{}=8`$–9 to $`Q^{}=13`$) is downhill in free energy.

It is also interesting to compare the free energy in the fly-casting corridor (Fig. 6b) with the free energy of folding for sequence S in isolation. For this purpose, we performed a simulation of S in isolation at $`T=ϵ/2.268k`$, where the native population (29%) is comparable to the binding probability at $`T=ϵ/2.6k`$ (35%). From Fig. 7, it can be seen that the free energy for the isolated chain is markedly different from that for the fly-casting corridor. In particular, we see that the bound-state minimum of the interacting chain is broader than the native minimum of the isolated chain; the isolated chain must reach $`Q^{}=11`$ before the free energy starts to decrease.

We also performed Monte Carlo-based kinetic simulations (see Sec. II) of the sequences S and U, using the same temperatures as in the thermodynamic runs. For each sequence, we carried out a set of 500 simulations, each containing $`5\times 10^6`$ Monte Carlo sweeps. The simulations were started from configurations obtained by short runs at a higher temperature ($`T=ϵ/2k`$). The resulting configurations from these preparatory runs (the starting configurations for the kinetic runs) showed a wide variation in both location and conformation for the moving chain. The averages of $`I`$ and $`Q`$ over this ensemble were about $`12`$ and 4–5, respectively.

Figure 8 shows the evolution of this ensemble of 500 systems with Monte Carlo time in the ($`I,Q`$) plane (not $`Q^{}`$). This parametric plot, with Monte Carlo time as a parameter, has the advantage of being independent of the difference in time scale between the two different temperatures used. To reduce noise, each data point represents an average over $`10^4`$ Monte Carlo sweeps, thus giving 500 data points in total. In Fig. 8, we also indicate the equilibrium values of $`I`$ and $`Q`$, which were obtained by a separate, very long simulation. The equilibrium values are given by $`(I,Q)=(6.1,12.1)`$ for S and $`(I,Q)=(4.9,11.1)`$ for U. Note that, as expected, S has a higher value of $`Q`$ and U has a lower value of $`I`$.

For sequence S, we see that the relaxation towards the equilibrium point has a clear two-step character, where the first step corresponds to folding (increasing $`Q`$) and the second step to binding (decreasing $`I`$). This behavior is indeed what one expects if docking is the preferred binding mode. For sequence U, there is no such clear distinction between folding and binding; $`I`$ and $`Q`$ evolve in a more correlated manner. Hence, the results of the kinetic simulations support the conclusions from the free-energy analysis.

## IV Summary

We have studied different binding mechanisms of polypeptide chains using the minimal HP model. The aim of this study is to understand the nature of an important aspect of protein-protein recognition and it can serve as a guideline for further, detailed studies of this complex phenomenon.

For small chain lengths, we found that the fraction of sequences that can bind to a given (fixed) target is comparable to the fraction of sequences that are stable in isolation. The overlap between these two sets of sequences is small; most binding sequences are unstable in isolation and fold upon binding to the target.

We then compared the binding behaviors of two related $`N=25`$ sequences. One of the sequences, S, had been optimized for high stability in isolation,Irback:02  and was found to have rigid docking as its preferred binding mode. The other, slightly mutated sequence, U, was found to have a bound structure identical to that of S, but to be unstructured in isolation. In these respects, sequence U is reminiscent of the engineered protein ZSPA-1 mentioned in the introduction. It turned out that sequence U tends to first attach to the target and then fold, which supports the conclusionShoemaker:00  that unstructured chains prefer binding through a fly-casting mechanism. The free-energy profile associated with this binding mode was found to lack high barriers and to exhibit a broad bound-state minimum, compared with the native minimum of sequence S in isolation. The observed difference in shape between these free energies suggests that by attaching itself to the target, the chain becomes able to fold more efficiently than it does in isolation. The free energy of docking obtained for sequence S contains, as expected, a significant entropic barrier, which the chain has to overcome in order to reach the bound state.

These free-energy profiles were obtained using two particular sequences in this model. Therefore, it should be stressed that the shape of the free energy was not considered at all when selecting these two sequences and the target structure; the goal was just to have two binding sequences with the same bound structure but different stability properties in isolation. As a result, we believe that the trends seen have some generality. Nevertheless, it is clear that it would be highly desirable to extend these calculations to other and longer sequences and to more realistic off-lattice models. Furthermore, it would be interesting to study the kinetics over a larger time interval; this requires a larger number of systems than 500 which we used, in order to keep the statistical errors under control.

### Acknowledgments

We thank Björn Samuelsson for valuable discussions and help with figures. This work was in part supported by the Swedish Research Council.


text:
# Phonons and structures of tetracene polymorphs at low temperature and high pressure

## I Introduction

Among molecular organic semiconductors, oligoacenes crystals represent a subject of increasing experimental and theoretical interest, because their high carrier transport properties make them likely candidates for applications in electronic and opto-electronic devices. Gundlach02 ; Butko03 ; Boer03  New techniques have been exploited for the growth of acene ultrapure single crystals Butko03 ; Boer03 ; Laudise98  or thin solid films, Voigt03  with the aim of obtaining high quality, well ordered crystalline samples. In fact, the absence of structural defects is crucial for the achievement of optimal performances in charge carrier mobilities. Boer03

Much effort has been devoted to clarify the polymorphism of pentacene. Venuti02 ; Brillante02 ; Valle03 ; Mattheus03 ; Siegrist01  Starting from all the published x-ray structures for crystalline pentacene, Siegrist01 ; Mattheus01 ; Campbell62 ; Holmes99  we computed the structures of minimum potential energy, and obtained two local minima of the potential energy, i.e., two different “inherent structures” of mechanical equilibrium. Venuti02  This behavior indicated that there were at least two different single crystal polymorphs of pentacene. The calculations predicted significant differences between the corresponding Raman spectra of the lattice phonons, which we checked experimentally, confirming the existence of two polymorphs. Brillante02  The correct identity of the samples, initially assigned only by matching experimental and calculated spectra, was thus verified directly with x-ray diffraction measurements. Finally, we obtained theoretical information on the global stability of the minima by systematically sampling the potential surface of crystalline pentacene. Valle03  We found that the two polymorphs correspond to the two deepest minima. Further deep minima with layered structures, which might correspond to the thin film polymorphs found to grow on substrates, were also predicted.

The existence of high temperature (HT), Campbell62 ; Holmes99  low temperature (LT) and high pressure (HP) polymorphs Prikhotko66 ; Vaubel70 ; Turlet73 ; Kolendritskii79 ; Jankowiak79 ; Sondermann85 ; Kalinowski78 ; Rang01  of tetracene has been known and studied in the past, and has been reported recently in studies on electronic transport in tetracene single crystals. Boer03 ; Boer04  Phase transitions for this system seem indeed to occur under variable conditions, depending on sample preparation, history and cooling speed. Boer03 ; Jankowiak79 ; Sondermann85 ; Boer04  As the transformations are sluggish, the sample can show a large temperature range in which more than one structure is present, and important hysteresis effects can be observed. Generally, this results in lowered carrier mobilities and shattering of the crystal upon cooling. Altogether, however, not much is known about the characteristics of the transitions and the nature of the polymorphs involved. For instance, it is not clear yet how many phases are actually formed at low temperature, or whether a LT phase does correspond to the HP one.

In this paper we address the issue of polymorphism in tetracene with the methods already successfully used for pentacene polymorphs. First, we provide lattice phonon Raman spectra obtained by means of a microprobe technique for two different phases of solid tetracene as a function of both temperature and pressure. Interfacing optical microscopy to Raman spectroscopy allows for a detailed mapping of the physical features of each crystalline sample and probe the conditions under which more phases can be simultaneously present. Secondly, experimental data are compared with the results of quasi harmonic lattice dynamics Ludwig67 ; Valle96 ; Valle98  (QHLD) calculations performed by using either the available crystallographic data, Campbell62 ; Holmes99 ; Sondermann85  or the theoretical predictions of the most stable crystal structures for this system, which were obtained by performing a systematic sampling of the potential energy surface, Valle04  as already done for pentacene.

## II Experimental

High temperature crystalline tetracene (HT structure) can be obtained by sublimation under vacuum of the commercial product (Aldrich) as thin platelets of the length of a few mm and thickness of tenths of $`\mu `$m. Samples of microcrystalline powder obtained by a fast sublimation process may instead contain domains of the low temperature (LT) structure as physical impurity. Sublimation in an inert atmosphere at reduced pressures (10–20 kPa of nitrogen or argon) at 493 K, with a procedure similar to that used to obtain polymorph II of pentacene, Brillante02  yielded directly larger amount of physically pure LT phase.

Raman scattering has been detected by using several laser lines, eventually selecting low energy excitation from a krypton laser tuned at 752.5 nm to minimize sample fluorescence, due perhaps to some residual chemical impurities. Boer03  Raman spectra above 1 GPa were anyway overlapped by strong fluorescence in all conditions. The spectra have been collected and analysed by the Jobin Yvon T64000 spectrometer equipped with a liquid nitrogen cooled CCD detector. Low temperatures $`T`$ down to 80 K were achieved in a conventional cryostat (Linkam HFS 91) with a temperature gradient of 10 K/min. High pressures $`p`$ up to 6 GPa were obtained in a LOTO diamond anvil cell, Loto  using perfluorocarbon as pressure medium. Pressures were measured with the ruby luminescence method. Piermarini75

Low $`T`$ and high $`p`$ cells were placed on the stage of the microscope (Olympus BX40) directly interfaced to the spectrometer. The use of 20x and 50x magnification objectives allowed for a spatial resolution of 2.2 and 1.1 $`\mu `$m, respectively, yielding the possibility to spatially check the physical purity of crystal polymorphs by mapping the lattice phonon profiles along the sample surface. Brillante02

## III Computational methods

The crystal structures and vibrational frequencies of tetracene have been calculated with the same procedure used for pentacene, Brillante02 ; Valle03 ; Masino02  following a well assessed treatment. Valle01 ; Girlando00  We first compute ab-initio molecular geometry, atomic charges, vibrational frequencies and cartesian eigenvectors of the normal modes for the isolated tetracene molecule. This is done with the Gaussian98 program Frisch98  (Rev. A.5), using the 6-31G(d) basis set combined with the B3LYP exchange correlation functional. Frisch98 ; Lee88  The vibrational frequencies are scaled by the factor of 0.9613 recommended Rauhut95 ; Scott96  for the combination of B3LYP and 6-31G(d).

The crystal total potential energy $`\mathrm{\Phi }`$ is given in terms of an atom-atom Buckingham model, Califano81  with Williams parameter set IV, Williams67  combined with an electrostatic contribution represented by a set of ab-initio atomic charges. We have chosen the potential derived charges, Frisch98  which describe directly the electrostatic potential.

The effects of temperature and pressure are accounted for by computing the structures of minimum Gibbs energy $`G(p,T)`$ with a QHLD method.Ludwig67 ; Valle96 ; Valle98  In this method, where the vibrational Gibbs energy of the phonons is estimated in the harmonic approximation, the Gibbs energy of the system is $`G(p,T)=\mathrm{\Phi }+pV+_ih\nu _i/2+k_BT_i\mathrm{ln}\left[1\mathrm{exp}\left(h\nu _i/k_BT\right)\right]`$. Here $`V`$ is the molar volume, $`_ih\nu _i/2`$ is the zero-point energy, and the last term is the entropic contribution. The sums are extended to all phonon frequencies $`\nu _i`$. Like pentacene, also tetracene exhibits ab-initio vibrational frequencies in the energy range of the lattice modes, and the coupling between lattice and intramolecular vibrations cannot be neglected. Valle01 ; Filippini84  To account for it, we adopt an exciton-like model, Girlando00 ; Califano81  where the interaction between different molecular coordinates is mediated by the intermolecular potential which depends directly on the atomic displacements. Since these correspond to the cartesian eigenvectors of the normal modes of the isolated molecule, we use the ab-initio eigenvectors and the scaled ab-initio frequencies. Intramolecular modes above 300 cm-1 are not taken into account, as the coupling is expected to be important only for low frequency modes. Masino02

## IV Temperature dependence of Raman spectra at constant pressure

For tetracene two distinct bulk crystal structures have been detected by x-ray diffraction experiments. At room $`p`$,$`T`$ conditions tetracene crystal was found to be triclinic, Campbell62  space group $`P\overline{1}`$ ($`C_i^1`$). The unit cell contains two independent molecules located on the (0,0,0) and ($`\frac{1}{2}`$,$`\frac{1}{2}`$,0) inversion sites. Thus, the factor group analysis for the lattice phonons at $`k=0`$ predicts six Raman active modes of $`A_g`$ symmetry and three IR active modes of $`A_u`$ symmetry. More recently Holmes et al. Holmes99  reported full data for crystalline tetracene at 183 K. As seen in section III, and unlike what found for pentacene, Venuti02 ; Brillante02  the two experimental structures of refs. Campbell62, ; Holmes99,  belong to a single polymorph. Finally, in 1985 Sondermann et al. Sondermann85  identified by x-ray diffraction a second triclinic polymorph, for which unit cell parameters at 140 K were given, but it was not possible to determine whether the space group was $`P\overline{1}`$ or $`P1`$.

Room $`T`$ Raman spectra of crystalline tetracene were reported in the late seventies by Jankowiak et al., Jankowiak79  who also reported Raman spectra as a function of $`T`$. Discontinuous changes in the temperature dependence of Raman lattice modes and the appearance of new phonon lines at 182 and 144 K were interpreted with the occurrence of two phase transitions, Jankowiak79  although some doubts could be harbored about the evidence of the high temperature (182 K) one. Also, a number of spectroscopic methods were used to test the occurrence of low temperature polymorphs in tetracene, Prikhotko66 ; Vaubel70 ; Turlet73 ; Kolendritskii79 ; Kalinowski81  yielding a quite large range of transition temperatures and conditions. The shattering of the sample upon cooling and large hysteresis effects upon heating the low $`T`$ structure have been observed as a consequence of the phase transition. Jankowiak79 ; Sondermann85

With he aim of rationalizing the situation, we started the experiments by measuring Raman spectra at room $`p`$, $`T`$ in the wavenumber range 20–300 cm-1 for samples grown by sublimation in a variety of ways, as described in section II. Depending either on the method of preparation or history of the samples, two different phonon patterns can be observed even at ambient conditions. To clarify the issue, we report in Figure 1 (bottom trace) the Raman spectrum at 298 K of the thin tetracene platelets typically obtained by sublimation under vacuum. In the same Figure (top trace) we also show the spectrum at the same temperature of samples grown in an inert atmosphere at reduced pressure. Whereas the former spectrum agrees with that reported Jankowiak79  by Jankowiak et al. for the HT structure (hereafter called polymorph I), the latter is instead found to overlap the spectra of samples recovered at ambient conditions after both LT and HP cycles, as it will be shown in the following. Therefore, these samples correspond to another phase, which will be called polymorph II. Most microcrystalline specimens (Figure 1, centre trace) display the typical bands of polymorph I with additional features of variable intensity, which can be attributed to different amount of polymorph II present as physical impurity.

The phonon wavenumbers for polymorphs I and II at room $`T`$ are given in Table 1. Since the number of Raman bands observed for polymorph II is always six, as for polymorph I, it is likely that also polymorph II belongs to the $`P\overline{1}`$ space group symmetry. As already revealed by the spectra of Figure 1, the phonon pattern of the two polymorphs differs especially in the lowest frequency region. Polymorph I displays three closely grouped bands in the range 42–58 cm-1, while polymorph II shows three evenly spaced bands in the range 38–73 cm-1. We remark that this finding closely resembles what already reported for polymorphs C and H of pentacene, Brillante02  respectively.

For the low temperature measurements we chose the thin platelets of physically pure polymorph I, to assure that the starting material for the temperature cycling was physically homogeneous and belonging to the structure thermodynamically stable at ambient conditions. Selected spectra recorded on decreasing temperature in the range 298–80 K are shown in Figure 2. No discontinuities were seen in our samples down to 140 K. At 130 K the abrupt appearance of the pattern typical of polymorph II is observed in the low frequency phonon region. The spectral changes are either accompanied or preceded by a cracking of the crystal. However, unlike what reported by Jankowiak et al., Jankowiak79  there is no hint of an intermediate crystal modification occurring in the range 180–140 K, even after repeated temperature cycling on several different specimens. The phase transformation is clearly completed at 80 K, where no features of the high temperature spectrum remain. The phonon wavenumbers of the 80 K spectrum of Figure 2 are reported in Table 1, and should be compared with the (incomplete) spectrum at 77 K of ref. Jankowiak79, .

Interesting information can be obtained by repeated temperature cycling, as shown in Figure 3, where we display a sequence of spectra of a single sample. After the temperature transition has occurred (Figure 3 a,b), a large hysteresis is documented by the persistence of polymorph II upon heating up to 298 K, where all spectral features of this phase are still retained on the time scale of the experiment, as shown in Figure 3 c. Note that the latter spectrum overlaps the top one of Figure 1. This hysteresis effect allowed us to measure the temperature dependence of the phonon bands for polymorph II by cooling the sample over the same $`T`$ range of polymorph I (Figure 3 d). Only by heating again polymorph II up to 320 K the conversion to polymorph I begins, although the process could be completed only by annealing at 400 K (Figure 3 e). The data will be compared with calculations in Section VI.

To summarize, two polymorphs have been clearly identified by the Raman analysis as function of temperature. Polymorph I is the most frequently grown, and it is the form stable at room $`T`$. Polymorph II is the form obtained by lowering temperature below 140 K. However, it can be obtained as a (metastable) phase also at room $`T`$, either by sublimation at 493 K at reduced pressure in an inert atmosphere or by bringing back to room $`T`$ samples cooled down to 80 K.

## V Raman spectra under pressure

High pressure induced transformations in tetracene crystals were observed as discontinuous changes of the Davydov splitting in the electronic absorption spectrum, Kalinowski78  spatial anisotropy of the magnetic field effect on fluorescence Kalinowski76  and recently by studying the photoconductivity of tetracene single crystals. Rang01  Mechanical stress induced by sample grinding was reported to produce a mixture of different phases. Jankowiak79 ; Sondermann85  So far, no Raman spectra were reported for the pressure induced transformation, and no satisfactory characterization was performed for the HP phase. Starting from polymorph I, we have tried to record Raman spectra as a function of $`p`$ over the range 0–6 GPa. However, in all samples and for all measuring conditions a strong fluorescence emission hides the Raman scattering above 1.4 GPa, and even in lower pressure regimes the detection of all bands of the spectrum turns out to be quite difficult, as shown in Figure 4. Raman spectra of our samples show that even at 1.0 GPa there is no sign of phase change. Therefore, the onset of the pressure-induced phase transition is well above the value (around 0.3 GPa) previously reported Kalinowski78 ; Rang01 ; Kalinowski76  with different experimental techniques. In any case, the investigation of the samples recovered at ambient $`p`$ clearly shows that a complete transformation has taken place, as the spectrum displays all the bands polymorph II and none of polymorph I (Figure 4). Therefore, polymorph II can be obtained from polymorph I also by pressure cycling, in addition to temperature cycling described in the previous Section.

## VI Computational results

We now describe the computational results, which help to definitely clarify the nature of the two tetracene polymorphs identified through Raman spectroscopy. In an early stage of our work, we considered separately the two complete experimental structures. Campbell62 ; Holmes99  Each structure was modeled starting from its experimental molecular arrangements, by replacing the experimental molecular geometries with the ab-initio one. We thus discovered that these two structures, measured at different $`T`$, actually map into the same potential minimum, having identical energies and unit cells (Table 2). Therefore they correspond to a single phase, stable at room $`T`$, which we identify with polymorph I. Accordingly, the subsequent calculations for the $`T`$ dependence of structure and dynamics for this polymorph were all performed starting from the data of Holmes et al. Holmes99  and ab-initio molecular geometry. In Table 2, we compare the lattice parameters of the experimental structures of refs. Campbell62, ; Holmes99, ; Sondermann85,  to the calculated parameters of the structures at the minimum of $`\mathrm{\Phi }`$ and at the minima of $`G(p,T)`$. The latter were calculated at ambient pressure, and at the same temperatures of the experiments. As the various structures of the literature were not reported in the same standard crystallographic frame, we directly compare the equivalent reduced cells. Santoro70

In Table 2 we also report the minimum $`\mathrm{\Phi }`$ lattice parameters of the structures which are theoretically predicted to correspond to the two deepest minima on the potential energy surface. These have been identified by performing a systematic search for the minima of the potential energy hypersurface. Following the methods used for pentacene Valle03 , we used a Quasi Monte Carlo sampling scheme to generate several thousands of different initial structures. Starting from each structure, we then minimized the total potential energy by adjusting the cell axes angles, positions and orientations of the molecules. The technical details of the calculations, together with information on the overall distribution of minima, will be reported in a separate paper. Valle04  Here we observe that the two deepest minima present triclinic lattice, space group $`P\overline{1}`$ ($`C_i^1`$), and differ mainly for the orientation of the two molecules in the unit cell. As shown in Table 2, polymorph I of refs. Campbell62, ; Holmes99, ; Sondermann85,  maps very accurately onto the deepest minimum of the potential energy surface. The table also shows that the minimum $`\mathrm{\Phi }`$ structure for polymorph I reproduces reasonably well the experimental lattice parameters, with residual differences $``$ 3% for the unit cell axes and angles. The computed cell volume is $``$ 5% smaller than the experimental one at room temperatures. This discrepancy, which decreases upon cooling, is partly due to the thermal expansion, totally neglected in the minimum $`\mathrm{\Phi }`$ calculations. In fact, as it can be seen from the results of the minimum $`G(p,T)`$ calculations, including vibrational effects brings the calculated volumes within 2% of the experimental ones, or better, and reproduces correctly the thermal expansion of phase I. The experimental volume expands by 1.7% from 180 to 298 K, to be compared with a calculated expansion of 2.1%. Finally, very good is the agreement for the experimental sublimation heat (unspecified phase) $`\mathrm{\Delta }_{\mathrm{sub}}H=34.4\pm 1.2`$ kcal/mol, Dekruif80  which is to be compared with the Gibbs energy calculated at 0 K, $`G(0)35.1`$ kcal/mol.

A totally reliable comparison between our computed structures and the measurements for polymorph II is not currently feasible, since only the experimental cell parameters, and no atomic coordinates, are given in ref. Sondermann85, . An accurate comparison is also impossible for the theoretical structures predicted in the same paper on lattice energy considerations. Note, however, that the prediction of a crystal symmetry lowering from $`P\overline{1}`$ to $`P1`$ is not supported by our Raman measurements.

In our calculations, we have found that the minimum $`G(p,T)`$ structure computed at 140 K starting from the second deepest minimum of the potential energy surface Valle04  favorably compares with the experimental cell parameters Sondermann85  of polymorph II at the same $`T`$. The comparison is reported at the bottom of Table 2, and although the agreement with the experiment is with no doubt worse than for polymorph I, the association appears justified, and is well supported by the calculated phonon spectrum discussed below.

The experimental Raman frequencies recorded as a function of $`T`$ for polymorphs I and II are compared to the corresponding minimum $`G(p,T)`$ calculations in Figure 5. The data for polymorph I are compared to the results obtained starting from the structure by Holmes Holmes99  et al.. Instead, the data for polymorph II are compared to those obtained by starting from the second theoretical deepest minimum. Valle04  Lattice and intramolecular modes are both shown. In Table 1 we report the calculated Raman wavenumbers for the two structures at 298 K, along with their experimental values. For polymorph II we also give the values at 80 K.

For both polymorphs the first six Raman modes are almost fully intermolecular in character, with negligible intramolecular contributions. In fact, the lowest ab-initio vibration of $`g`$ symmetry in tetracene is calculated at 146.5 cm-1, and is weakly coupled only with the highest frequency lattice mode, which displays a intramolecular contribution of about 10% at 80 K. Therefore, in the Raman spectra tetracene behaves as an apparent rigid-body while the coupling is expected to be important for ir active modes of $`u`$ symmetry above 100 cm-1, as already remarked by Filippini and Gramaccioli. Filippini84

The temperature dependence of the phonon frequencies calculated for phases I and II (Figure 5), also agrees well with the corresponding experimental results. As expected, both experiments and calculations shows that varying the temperature affects the low frequency lattice modes much more than the purely intramolecular modes above 150 cm-1. The differences between the experimental and computed temperature dependence of the frequencies, especially noticeable for the high frequency modes, are attributed to the anharmonic frequency shifts,Valle83  neglected in these calculations, and to defects in the potential model.

To summarize the results of the phonon dynamics analysis, we point out again (cf. section IV) that the patterns of the experimental frequencies for the lattice modes of polymorphs I and II at room $`T`$ are clearly distinguishable, and are well matched by the patterns calculated for the experimental structure Holmes99  and for the second theoretical deepest minimum, respectively. In particular, the computations correctly predict the three closely spaced modes around 48 cm-1 for polymorph I, and reproduce the more widely spread bands between 38 and 73 cm-1 for polymorph II. Therefore we can associate polymorph I and II with the two deepest minima found in the potential energy surface.Valle04  Our calculations actually predict that polymorph I is the more stable (at 0 K) and denser phase, whereas experimentally the opposite is true. Moreover, we do not find theoretical evidence of the phase transition (crossing of $`G`$ values) as function of either $`T`$ or $`p`$. The same kind of problem has been found in calculations for pentacene.Masino02 ; Verlaak03  We could fine tune the atom-atom potential to yield the correct phase ordering. On the other hand, the calculated energy difference between the two phases is very small, less than 0.5 kcal/mole, which is the typical accuracy of this kind of calculations. We therefore think that attempts to improve the potential would be unjustified at the present stage.

## VII Discussion and conclusions

In this paper we have explored the $`p`$,$`T`$ phase diagram of crystalline tetracene. By combining Raman spectroscopy with computational methods, we have clarified several issues related to the crystalline phases of tetracene. Tetracene crystallizes into two different polymorphs, polymorph I and II, the former being the most frequently grown phase, stable at ambient conditions. Whereas the crystalline structure of polymorph I is well known, Campbell62 ; Holmes99  we suggest a likely structure for polymorph II, stable at low $`T`$ and high $`p`$. The structures of polymorph I and II are very similar to the structures of polymorph C and H of pentacene, respectively.Masino02  All structures are triclinic, space group $`P\overline{1}`$ ($`C_i^1`$), with two molecules per unit cell residing on symmetry unrelated inversion centers. In polymorph I and C, the long molecular axis is roughly pointing along the  direction ($`c`$ axis), in polymorphs II and H, which are the denser phases, the molecules are more inclined, and point towards the \[1,1,$`1`$\] direction. Raman spectroscopy in the lattice phonon region has indicated that the denser polymorph II of tetracene is stable only at low temperature (below 140 K) or high pressure (well above 1 GPa), whereas polymorph H of pentacene is the stable phase at ambient pressure. No additional phases have been detected, so early reports may have been affected by temperature induced strains and sample impurities/imperfections.Kalinowski78

Finally, we have verified the possibility of obtaining either or both tetracene polymorphs at ambient $`p`$,$`T`$ conditions, depending on sample preparation. Obtaining pure polymorph II at ambient conditions can be very important, as we expect that a denser phase should exhibit larger bandwidths and mobilities. The two tetracene phases are very similar in energy, and as it happens for pentacene, crystalline samples may show phase inhomogeneities. The two polymorphs can be easily identified through Raman spectroscopy in the lattice phonon region, and possible phase inhomogeneities are detectable this way. Raman spectroscopy thus represents a convenient and reliable tool for checking crystal quality, contribuiting to improve the performances of tetracene-based devices.

###### Acknowledgements.

Work supported by the Consorzio Interuniversitario Nazionale per la Scienza e Tecnologia dei Materiali (I.N.S.T.M., PRISMA2002 project) and by the Italian Ministero Istruzione, Università e Ricerca (M.I.U.R., FIRB project).


text:
# A Mechanism for Polar Protein Localization in Bacteria

## I Introduction

The polar localization of proteins is a vital component of the sophisticated internal organization of bacteria, and is important for many processes ranging from chemotaxis to accurate cell division Shapiro ; Errington ; Lutkenhaus ; Shapiro2 . However, how such polar localization is achieved has remained largely unknown. The purpose of this report is to analyze one possible mechanism for localization in a well–studied experimental context: the MinCD/DivIVA system regulating cell division positioning in the Gram–positive bacterium B. subtilis Cha ; Edwards1 ; Marston1 ; Marston2 ; Karoui ; Muchova ; Harry ; Hamoen . We will then discuss how these ideas might be applied more generally to the generation of polar localization/asymmetry in bacteria.

The dynamics of bacterial cell division have been intensively studied, particularly in B. subtilis Cha ; Edwards1 ; Marston1 ; Marston2 ; Karoui ; Muchova ; Harry ; Hamoen  and E. coli Raskin1 ; Raskin2 ; Hu1 ; Yu ; Rowland ; Fu ; Rothfield ; Hale ; Hu2 ; Suefuji ; Shih1 , where the polar localization of proteins is known to be crucial. In both cases, the location of the division site is determined by the tubulin–like GTPase FtsZ which assembles into a ring–shaped structure just under the cell membrane at midcell. In E. coli accurate targeting of the FtsZ ring to midcell is achieved both by the inhibitory effect of the nucleoids (“nucleoid occlusion”) and through the MinCDE proteins Yu . MinC and MinD form a complex on the cell membrane, which undergoes coherent pole–to–pole oscillations, with a period of around a minute, driven by MinE. As a result of the oscillations, the time–averaged MinCD concentration is minimized at midcell and maximized at the cell poles. Since MinC inhibits FtsZ polymerization, the FtsZ ring assembles at midcell where the average MinC concentration is lowest. In this way, through the temporary sequestration of MinCD at the cell poles, the midcell can be accurately targeted for cell division.

Recently, several groups have employed simple reaction–diffusion mathematical models (either deterministic Howard1 ; Meinhardt ; Kruse  or stochastic Howard2 ) to try to understand the MinCDE dynamics (see also Refs. Meinhardt2 ; Hunding  for related models). The principal result of this modeling effort is that the oscillations are a spontaneous self–organized phenomena arising from an intrinsic instability of the reaction–diffusion dynamics. A vital element of this instability is the varying diffusion constants of the MinCDE proteins: high in the cytoplasm; low when polymerically bound to the membrane. Crucially, the oscillations form spontaneously as a result of this disparity of the diffusion constants and of the nonlinear binding/unbinding interactions between the proteins. As a result, the oscillatory dynamics does not require any pre–existing topological markers to distinguish various locations within the cell, as had previously been supposed.

Very recently, improved fluorescent imaging has revealed that the E. coli Min proteins polymerize into stable helices on the cytoplasmic membrane Shih1 . This phenomenon is not explicitly included in any of the existing mathematical models, which may therefore require some amendment. Nevertheless, the models generally agree well with fluorescent imaging experiments, especially in the observation of periodic MinCDE/FtsZ bands in filamentous mutants. These results are very difficult to understand other than as a pattern–forming reaction–diffusion instability. Hence, we believe the essentials of a reaction–diffusion mechanism underpinning the MinCDE oscillations are likely to be correct. Although the formation of the Min proteins into filaments has only been observed in E. coli (and also in vitro Suefuji ), this observation may have important implications in other bacteria. Already in B. subtilis the shape–determining proteins MreB and Mbl are known to form helical cables that could comprise cytoskeletonal–like elements Errington ; Jones .

Despite possessing two of the three same proteins (MinC and MinD, but not MinE), B. subtilis appears to use a quite different mechanism to achieve central placement of the FtsZ ring in vegetative cells Cha ; Edwards1 ; Marston1 ; Marston2 ; Karoui ; Muchova ; Harry ; Hamoen . In B. subtilis, MinCD do not oscillate, but instead are localized to both cell poles, where they are anchored by the protein DivIVA (which is structurally unrelated to MinE). Also involved are the (initially) homogeneously distributed proteins ZapA, which promotes Gueiros , and EzrA, which inhibits Levin1 , FtsZ ring assembly. Quite how DivIVA is able to recognize the cell poles has remained unresolved and forms the main object of this report. With MinCD localized to the cell poles, the FtsZ ring is again unable to nucleate there. However, it is likely that additional mechanisms are used to accurately locate the FtsZ ring to midcell Levin2 , since the MinCD/DivIVA system functions primarily to block polar division Migocki . These additional (and so far uncharacterized) mechanisms may be especially important as nucleoid occlusion is ineffective as a means of positioning in B. subtilis Regamey . Despite this, the precision of midcell division is the same in E. coli and B. subtilis Migocki , even though their regulatory mechanisms differ considerably.

Polar protein localization is not only restricted to positioning of the division site. For example in Caulobacter crescentus, polar asymmetry is exploited to ensure that daughter cells are dissimilar: one a swarmer cell with a polar flagellum, and one a stalked cell with a polar stalk that can adhere to a source of food Shapiro ; Shapiro2 . In each case, the development of a cell with a polar stalk or flagellum is critically dependent on spatiotemporally controlled protein localization, and on the subsequent proteolysis and release of these components at the appropriate time. Other examples of polar localization include the polar clustering of chemoreceptors, in both E. coli Maddock ; Lybarger ; Maki  and C. crescentus Alley , and also the polar localization of the virulence protein IcsA Robbins , which promotes undirectional actin–based motility of the bacterium Shigella flexneri.

## II Polar protein localization in B. subtilis

Polar division in vegetative cells of B. subtilis is blocked by the presence of membrane localized MinC. MinC is normally a cytoplasmically localized protein, however, in the presence of MinD the two interact strongly and bind to the membrane Marston2 . Polar localization is facilitated by the action of the protein DivIVA, which acts to anchor MinCD to the poles Marston1 ; Marston2 . Targeting of MinCD to the poles, is dependent on DivIVA, as in filamentous cells lacking DivIVA, MinCD was partly distributed diffusively along the cell Marston1 ; Marston2 . Furthermore the localization and activity of MinD is believed to require ATP binding and possibly hydrolysis Karoui . However, MinCD can also (MinC more strongly than MinD) target past/present division locations independently of DivIVA Marston1 ; Marston2 . Similarly in a minD mutant, DivIVA simply localizes to division sites regardless of location, again demonstrating an affinity for some component of the cell division apparatus Marston1 . Thus in wild–type cells, MinCD/DivIVA localize to the poles (old division sites), blocking cell division there. However, before septal constriction occurs and the cell divides, MinCD/DivIVA assemble approximately simultaneously at midcell Marston1 , and are retained there after division Marston1 ; Marston2 . In this way MinCD is delivered to the new poles of the two daughter cells, meaning that polar division is immediately prevented at the new poles.

This observation immediately provides one possible explanation for polar localization: since the cell division apparatus marks the location of the new cell poles in the daughter cells, if the relevant proteins can localize to a component of the division apparatus, then the cell pole will continue to be marked in the daughter cells. However, one may question whether this mechanism is really sufficiently robust: if MinCD/DivIVA fail to correctly localize to the site of division, then the poles of the resulting daughter cells will remain unmarked, meaning that division in those cells would be corrupted. It would certainly therefore be desirable to have other means of identifying the cell poles. Experiments in outgrowing spores of B. subtilis indicate that this is indeed the case Harry ; Hamoen . MinCD/DivIVA are targeted to the cell poles very rapidly after spore germination, even in cells lacking FtsZ. Thus, polar localization can occur even in germinating spores which do not possess pre–existing division apparatus left over from prior division events.

So the question arises, what other mechanism could be responsible for the polar localization? One possible answer is that a self–organized reaction–diffusion instability exists, similar to that already studied in E. coli. Of course, for B. subtilis, the instability would have to be to a steady–state distribution where the MinCD/DivIVA do not oscillate (i.e. in linear stability analysis one would expect a real rather than a complex growth eigenvalue). With appropriate interactions and diffusion constants, such an instability could certainly exist. However, one of the key predictions of such an instability is the existence of a characteristic wavelength, depending on the microscopic parameters of the model, but independent of the geometry Howard1 ; Meinhardt ; Kruse . Thus in long filamentous cells of E. coli, several wavelengths can be fit into a single cell and one therefore expects to see multiple oscillating bands. The observation of this phenomena in filamentous E. coli Raskin1  is powerful evidence in favor of the reaction–diffusion mechanism. The situation in B. subtilis is, however, rather different: in filamentous cells lacking FtsZ, MinCD/DivIVA are primarily localized to the cell poles Marston1 ; Marston2 , without any regularly spaced intermediate bands in between, as one would expect for a reaction–diffusion instability. One could argue that the affinity of MinCD/DivIVA for the pre–existing (polar) division apparatus from previous division events might mask the existence of any intermediate bands. However, in outgrowing spores, which have no pre–existing division apparatus, there are still no intermediate bands, even in long filamentous mutants (up to about $`13\mu m`$ long) Harry ; Hamoen . Thus we see no evidence for the existence of a characteristic wavelength. Of course, it is just possible that the characteristic wavelength is simply very long, too long to have been seen in past experiments; but, there is currently no hint of this.

Instead we propose that it is a geometric effect in combination with reaction–diffusion dynamics which marks the end of the cell. The essence of our approach is to assume that MinCD has a reduced rate of membrane binding at the cell poles. If, for example, MinCD polymerize into filaments, as in E. coli Shih1 , this might be due to the increased curvature of the cell poles interfering with membrane polymerization. Since DivIVA may then bind to the edges of MinCD clusters and so stabilize them, the MinCD (and DivIVA) concentrations end up being maximized close to the cell poles, even though MinCD initially has a reduced rate of binding there. Similarly, when the cell divides a region of high curvature is created at midcell, again attracting the MinCD/DivIVA proteins independent of the division apparatus. In the next section, we define this model more precisely and report on the results of numerical simulations.

## III The model

Previous modeling work on the MinCDE system has demonstrated the effectiveness of simple deterministic reaction–diffusion equations Howard1 ; Meinhardt ; Kruse . Accordingly we follow a similar modeling approach here. The dynamics of MinC is similar to that of MinD, so, for simplicity, we drop MinC from our model. The fact that MinD can localize normally in minC mutant cells also supports this approach Marston2 . Initially, we also assume a “conservative” model, where the MinD/DivIVA are able to bind/unbind from the membrane but are not degraded in the process. This is again similar to some of the models used for E. coli Howard1 ; Kruse ; Howard2 . The relevant set of one–dimensional reaction–diffusion equations are:

$`{\displaystyle \frac{\rho _D}{t}}=D_D{\displaystyle \frac{^2\rho _D}{x^2}}{\displaystyle \frac{\lambda _D(x)\rho _D}{(1+\lambda _D^2\rho _d^2)}}+{\displaystyle \frac{\lambda _d\rho _d}{(1+\lambda _d^{}\rho _a)}}`$
$`{\displaystyle \frac{\rho _d}{t}}=D_d{\displaystyle \frac{^2\rho _d}{x^2}}+{\displaystyle \frac{\lambda _D(x)\rho _D}{(1+\lambda _D^2\rho _d^2)}}{\displaystyle \frac{\lambda _d\rho _d}{(1+\lambda _d^{}\rho _a)}}`$
$`{\displaystyle \frac{\rho _A}{t}}=D_A{\displaystyle \frac{^2\rho _A}{x^2}}{\displaystyle \frac{\lambda _A\rho _A\left|\frac{\rho _d}{x}\right|}{(1+\lambda _A^2\rho _a^2)}}+{\displaystyle \frac{\lambda _a\rho _a}{(1+\lambda _a^2\rho _a^2)}}`$
$`{\displaystyle \frac{\rho _a}{t}}=D_a{\displaystyle \frac{^2\rho _a}{x^2}}+{\displaystyle \frac{\lambda _A\rho _A\left|\frac{\rho _d}{x}\right|}{(1+\lambda _A^2\rho _a^2)}}{\displaystyle \frac{\lambda _a\rho _a}{(1+\lambda _a^2\rho _a^2)}}.`$
Here, the first equation is for the cytoplasmic MinD density $`\rho _D`$, the second for the membrane bound MinD density $`\rho _d`$, the third for the cytoplasmic DivIVA density $`\rho _A`$, and the fourth for the membrane bound DivIVA density $`\rho _a`$. In the above model, we assume that MinD binds to the membrane with rate $`\lambda _D`$, in agreement with green fluorescent protein (GFP) fusions which indicate that MinD is a membrane–localized protein Marston1 ; Marston2 . However, one would also expect that the rate of binding is reduced as the MinD membrane density increases. This accounts for the $`\lambda _D^{}`$ term. We also assume that MinD can unbind spontaneously ($`\lambda _d`$ term), but that this is suppressed by DivIVA ($`\lambda _d^{}`$ term) Marston1 ; Marston2 . It is also possible to insert a factor $`(1+\lambda _d^{\prime \prime }\rho _d)^1`$ in this term, which would model increased stabilization of membrane MinD from clustering. However, we have found that our results were not significantly affected by this term, which we have therefore left out for simplicity. Notice that the DivIVA functions to stabilize the membrane MinCD, hence its role is very different from the E. coli MinE protein which destabilizes membrane MinCD. This difference is partly responsible for the structure of the equations being rather different in this case from the earlier models for E. coli MinCDE Howard1 ; Meinhardt ; Kruse . DivIVA is also believed to be an intrinsically cytoplasmic protein Cha . Hence, in order to bind to the membrane in outgrowing spores Harry , (without pre–existing division apparatus), it presumably requires the presence of MinD. We assume here, similar to Ref. Karoui , that it binds to (and then stabilizes) the edges of MinD clusters. This could correspond to DivIVA “capping” the ends of MinCD filaments Karoui . This process would then account for the $`\lambda _A`$ gradient term $`|\rho _d/x|`$, which is maximized at the edges of a cluster. Note that the gradient term has a modulus sign, as the rate of binding should not depend on the sign of the gradient. Again we assume that the DivIVA binding is suppressed as the membrane DivIVA density increases ($`\lambda _A^{}`$ term). Finally DivIVA can spontaneously unbind ($`\lambda _a`$ term), but this is suppressed by DivIVA itself ($`\lambda _a^{}`$ term). This indicates that DivIVA itself oligomerizes/forms clusters for maximum stability, for which there is also experimental evidence Muchova .

The above assumptions are fairly minimal and are all closely based on the available experimental evidence of GFP fusions and deletion mutants Cha ; Edwards1 ; Marston1 ; Marston2 ; Karoui ; Muchova ; Harry ; Hamoen . As we have said, we then add one extra ingredient, namely a position dependent membrane binding $`\lambda _D(x)`$ for the MinD, where $`\lambda _D(x)`$ is reduced in regions close to the cell ends (see below).

For our computer simulations, we assume hard wall boundary conditions (i.e. $`\rho _i/x=0,i=\{D,d,A,a\},`$ at the boundaries) and we numerically integrate the above equations, with lattice size $`dx=0.02\mu m`$ and time increment $`dt=1\times 10^5s`$. We start from random initial conditions, where the MinD and DivIVA are randomly distributed throughout the cytoplasm. We begin the simulation with a total of $`50`$ lattice sites (bacterial length $`1\mu m`$), with initial densities of $`1000\mu m^1`$ for MinD and $`750\mu m^1`$ for DivIVA. The actual concentration levels for B. subtilis have not been accurately measured, so we assume concentration levels similar to those observed in E. coli Shih2 . During the course of the simulation, we then allow the bacterial length to increase, modeling a growing bacterium, with growth rate of $`4\mu m/\mathrm{hour}`$. In other words, every $`1.8\times 10^5`$ time steps an extra site is inserted, randomly located. Note, however, that the extra site is never inserted at the inert cell ends Daniel  (defined here by regions where $`\lambda _D(x)`$ is reduced, see below). As the bacteria grows in the simulation we add extra MinD/DivIVA, randomly localized in the cytoplasm, to ensure that the average densities remain constant (i.e. the total number of MinD/DivIVA proteins is proportional to the cell length). There is currently no information on whether the average densities do actually remain constant as the cell grows, but the assumption used here is the simplest possible. Nevertheless, precise regulation of the densities turns out not to be particularly important. Considerable overproduction eventually leads to the MinCD/DivIVA proteins saturating the poles and spilling over into the reminder of the cell length. Underproduction simply translates to less MinCD/DivIVA being present along the cell length in between the poles. The important point is that the bacterium has to manufacture sufficient MinCD/DivIVA to ensure that the cell poles are always well populated with these proteins, but not so much that the proteins spill over to midcell leading to division inhibition. Note that in our simulations, we exclusively model the case of outgrowing spores without pre–existing division apparatus. Clearly, modeling protein localization to pre–existing binding sites in wild–type vegetative B. subtilis cells is straightforward, and not particularly revealing.

In our simulations we use the following parameter set: $`D_D=2.0\mu m^2s^1`$, $`D_d=0`$, $`D_A=2.0\mu m^2s^1`$, $`D_a=2\times 10^4\mu m^2s^1`$, $`\lambda _D=10s^1`$, $`\lambda _D^{}=0.001\mu m`$, $`\lambda _d=3.0s^1`$, $`\lambda _d^{}=0.2\mu m`$, $`\lambda _A=0.03\mu m^2s^1`$, $`\lambda _A^{}=8\times 10^4\mu m`$, $`\lambda _a=2.0s^1`$, $`\lambda _a^{}=0.005\mu m`$. Furthermore, for the closest $`10`$ sites to both ends of the bacteria, the binding rate $`\lambda _D`$ is reduced by a factor of $`2`$ (our results are robust to changes in this suppression factor). It is also possible to generalize the model so that the MinD/DivIVA proteins are degraded at a (slow) rate, whilst increasing the production rate to keep the average densities constant. We have verified that this modification also makes no qualitative difference to our results. The values for the above cytoplasmic diffusion constants are of the same order as that measured (for different proteins) in the cytoplasm of E. coli Elowitz . As for the membrane diffusion constants, we have used similarly small values to those used in E. coli Howard1 ; Howard2 . Previously this was an important component in the spontaneous reaction–diffusion instability; this is, however, not the case here: the small membrane diffusivities are now necessary to preserve the polar density concentration maxima and to prevent the concentrations from equalizing everywhere along the cell membrane. For the remaining parameters there is no experimental guidance, but we have taken care to check that our results are robust to variations in these parameters.

## IV Results: polar protein localization

In Figure 1, we show the result of our computer simulation of the spore outgrowth experiment. We see very similar results to those observed experimentally Harry ; Hamoen . Starting from random initial conditions, the MinD/DivIVA rapidly (within a few minutes) localize to the poles of the cell. This agrees with experiments indicating that MinD/DivIVA localize to the cell poles very early in spore germination Harry ; Hamoen . In our model the mechanism for this localization follows from the reduced rate of MinD binding to the cell ends. This reduced polar value for $`\lambda _D(x)`$ is absolutely critical for polar localization as it sets up a pronounced membrane MinD density gradient close to the cell poles. The DivIVA then binds to these positions of increased gradient, stabilizing the MinD/DivIVA. In this way the MinD/DivIVA concentrations end up being maximized close to the cell ends. In fact, the concentrations are not maximized precisely at the poles, where there is actually a dip (see Figures 2 and 3). However the polar concentrations still end up being much higher than the densities well away from the cell ends. As the cell continues to grow, the polar regions become saturated with MinD/DivIVA and the newly produced proteins cannot localize there. Instead the proteins localize to random membrane patches along the cell length. Furthermore, these patches appear to propagate along the cell length (see Figure 1). This movement results from the coupling of the DivIVA to the edges of MinD clusters, slightly away from the local MinD density maxima. Since MinD is stabilized by the DivIVA this encourages MinD to rebind slightly away from its current density maxima. This in turn shifts the position of the maximum MinD density gradient, thus also encouraging the DivIVA to rebind in a shifted location. In this way the coupled MinD/DivIVA density maxima both tend to propagate along the membrane, as shown in Figure 1. It would be interesting to see if such movement can be observed in experiments. A snapshot of these patterns can be seen in Figure 2, where we show the MinD/DivIVA membrane density profiles at length $`10\mu m`$. The DivIVA pattern is very similar to that observed experimentally in filamentous spore outgrowth experiments Harry ; Hamoen , with sharp polar peaks and less intense randomly distributed membrane patches along the cell length. Note that this disordered pattern is seen everywhere inside the simulated outgrowing spore (including close to the cell ends) if $`\lambda _D(x)`$ is a constant, emphasizing again that polar protein localization is critically dependent on a position–dependent MinD binding rate.

More quantitative data, with profiles of the membrane MinD/DivIVA densities at different times are shown in Figure 3(a,c,e). In Figure 3(b,d,f) we show examples of the membrane density averaged over $`90s`$ as a function of position for various total cell lengths. A $`90s`$ average is chosen as, at least in E. coli, this is the rough time over which the FtsZ averages in locating the FtsZ ring Raskin1 . As can be seen, the MinD concentration is nicely minimized at the cell center for shorter cell lengths. However, as the cell length increases, the precision with which the center can be identified is considerably reduced, with the averaged signal fluctuating considerably at $`8\mu m`$. This data is consistent with evidence that the MinCD/DivIVA system is primarily responsible for blocking polar division, rather than in accurately placing midcell division Migocki . Notice also how the MinD is less tightly localized at the cell poles; its signal decreases more slowly than for the DivIVA towards the center of the cell. This also agrees with experiment Marston1 .

## V Results: asymmetric protein localization

Up to this point, we have modeled symmetric protein localization, such as that seen in the MinCD/DivIVA system in B. subtilis. However, it is also interesting to ask whether one can modify the model to obtain asymmetric localization, i.e. to one pole only. This modification turns out to be quite straightforward: in essence one merely has to remove the $`\lambda _D^{}`$ and $`\lambda _A^{}`$ terms which restrict membrane binding when the membrane densities are already high. In this case, starting from random initial conditions, the dynamics rapidly localizes the proteins to both cell poles. However, one randomly chosen pole will, due to fluctuations in the random initial conditions, gain a slightly higher density than the other. Since higher densities yield greater stability for the membrane aggregations (due to the $`\lambda _d^{}`$ and $`\lambda _a^{}`$ terms), and since there are no longer saturation effects, the favored pole, with an initially slightly higher concentration, will continue to attract additional proteins from the other pole, via cytoplasmic diffusion and then polar membrane capture (cf. Ref. Rudner ). Of course this process will also occur in reverse from the favored to the less favored pole, but because of its higher density, the net rate of accumulation will be towards the pole with the higher initial concentration. In this way, the less favored pole is eventually drained of its protein population. Starting from randomly distributed cytoplasmic MinD and DivIVA, with average densities $`1000\mu m^1`$, and $`750\mu m^1`$, respectively, an appropriate parameter set is: $`D_D=2.0\mu m^2s^1`$, $`D_d=0`$, $`D_A=2.0\mu m^2s^1`$, $`D_a=2\times 10^4\mu m^2s^1`$, $`\lambda _D=10s^1`$, $`\lambda _D^{}=0`$, $`\lambda _d=3.0s^1`$, $`\lambda _d^{}=0.001\mu m`$, $`\lambda _A=0.012\mu m^2s^1`$, $`\lambda _A^{}=0`$, $`\lambda _a=1.0s^1`$, $`\lambda _a^{}=1.2\times 10^4\mu m`$. In this case localization to a (randomly selected) pole occurs in about $`800s`$ (data not shown). Alternatively, replacing some of the above parameters with $`\lambda _d^{}=1\times 10^4\mu m`$, $`\lambda _A=0.002\mu m^2s^1`$, $`\lambda _a=1.6s^1`$, and $`\lambda _a^{}=1.6\times 10^4\mu m`$, unipolar localization occurs far more rapidly, within about $`40s`$. A process analogous to this may be happening in asymmetric division in sporulating B. subtilis Ben-Yehuda1  (see below), although in that case the protein SpoIIE may also play a crucial role Lucet .

## VI Discussion

In this paper, we have proposed a mechanism for the symmetric or asymmetric polar localization of proteins, relying on a combination of geometric information and reaction–diffusion dynamics. This mechanism could be particularly important for polar protein localization in cells growing out from spores, where there are no pre–existing division apparatus from prior division events. However, even in vegetative cells, where the poles could be marked by the division apparatus from prior divisions, the proposed mechanism could be useful for increased reliability (and redundancy) in locating the poles. We have focused particular attention on the MinCD/DivIVA system in B. subtilis, where our results agree well with experiments on outgrowing spores. We have also made a number of testable assumptions: very low diffusion constants for membrane–bound MinCD/DivIVA; DivIVA binds to the edges of MinD clusters; and membrane binding of MinD is heavily influenced by geometric effects (possibly by the high curvature of the cell poles restricting MinD membrane polymerization). We should point out, however, that experiments in round cell mutants have indicated that simple curvature of the cell membrane is insufficient to elicit DivIVA targeting Edwards2 . Nevertheless, this does not contradict our results: according to our model, what is important is that there be a difference in curvature: one part of the cell with low curvature, another part (the cell poles) with high curvature.

As we mentioned in the introduction there are several other regulatory systems in bacteria with similar polar localization properties. For example, the MCP chemotaxis receptor proteins (in E. coli) Maddock ; Lybarger ; Maki  and McpA (in C. crescentus) Alley  are both polar localized. Furthermore, filamentous mutants in both cases show an intriguing pattern of localization: the proteins are still polar localized, but in addition there is a seemingly randomly distributed component along the filament length between the two poles. This pattern is qualitatively reminiscent of the pattern seen (and modeled in this paper) for MinCD/DivIVA in filamentous B. subtilis mutants grown from spores Harry ; Hamoen . Although, for the chemotaxis cases, it may simply be that the proteins have saturated polar binding sites and are then “spilling” over into the remainder of the cell. However it is again possible that there may be a geometric component to the localization, similar to that proposed in this report. This may be even more likely in the case of chemotactic proteins in outgrowing spores. Already the chemotaxis protein McpB is known to be polar localized in vegetative B. subtilis Kirby . It would be very interesting to study its localization properties in outgrowing spores without pre–existing division apparatus. In addition, it has been shown that DivIVA localizes to division sites both in E. coli and, unexpectedly, in the fission yeast Schizosaccharomyces pombe Edwards2 . However, this may again be due to binding to division apparatus proteins (possibly actin in S. pombe).

Given that the geometric properties (high curvature) of the cell poles may be important in influencing MinD binding, it is, of course, important to understand how the cell shape is maintained. In B. subtilis, cell elongation and maintenance of the cell cylinder is achieved by Mbl–directed dispersed helical insertion of cell wall material throughout the cylinder, while the hemispherical cell poles are effectively inert Daniel . In vegetative cells, where a midcell septum is formed, it is plausible to see how the hemispherical end is constructed and thus how overall cell shape maintained. However, in spore outgrowth, this is not so obvious. B. subtilis spores themselves are already ellipsoidally shaped, so the initial state for spore germination already has the ends marked via high curvature. However quite how this shape is produced in the earlier sporulation event is rather unclear.

It is also intriguing to note that, although MinC and MinD are used by both B. subtilis and E. coli, the two bacteria appear to use entirely different physical mechanisms to exclude division from polar sites. E. coli uses a spontaneous reaction–diffusion instability, whereas geometric effects (reduced polar binding of MinD) and associated reaction–diffusion kinetics, appear to be employed in B. subtilis. Interestingly, it is actually a reduced rate of MinD membrane binding at the cell poles which, in our model, eventually leads to polar localization due to DivIVA binding/stabilizing the edges of MinD membrane clusters. It would be interesting to further study these mechanisms with their evolutionary origin in mind, especially since both employ reaction-diffusion dynamics. Although we believe that the mechanism outlined in this paper provides a very plausible scenario for polar protein localization, other mechanisms are still possible. For example, one could imagine that the DivIVA makes the poles more attractive for MinCD polymerization, perhaps by altering the structure of the MinCD filaments so that they prefer a higher curvature region (see also Ref. Harry ). Given that the model described in this paper provides several testable predictions, experiments should be readily able to discriminate between these possible models.

One may also speculate on why E. coli and B. subtilis have differing mechanisms. One possible explanation is that B. subtilis must have extra flexibility in its cell division regulation as, unlike E. coli, it sporulates when confronted with harsh environmental conditions. In the process of forming a spore, the division septum is relocated from midcell to one of the cell poles. Hence B. subtilis must have the ability to undo its polar division block. All that would be required is for another protein to bind to DivIVA, displacing the MinCD, thus freeing up the polar sites, as suggested in Ref. Ben-Yehuda2 . A suitable candidate for such a displacing protein has recently been shown to be RacA Ben-Yehuda2 . The MinCD presumably then relocates to the cell center, forcing the FtsZ ring to spiral outwards to the cell poles, as observed in Ref. Ben-Yehuda1 . As only a single pole is then selected for spore formation, the cytokinetic proteins should localize asymmetrically, possibly in a similar way to that discussed above.

Finally, we mention that the repertoire of division site positioning is likely still not exhausted by the MinCDE and DivIVA systems analyzed here and in previous modeling work Howard1 ; Meinhardt ; Kruse ; Howard2 . C. crescentus is able to divide asymmetrically with high precision even though it lacks all of the MinCDE and DivIVA proteins Shapiro2 . How this is achieved remains unknown and hints that there may be yet other mechanisms capable of providing accurate positional information and hence of regulating subcellular architecture.

Acknowledgements

This work was supported by The Royal Society.