A simple mechanism for controlling the onset and arrest of collective oscillations in genetic circuits
Abstract
We study a system of dynamical units, each of which shows excitable or oscillatory behavior, depending on the choice of parameters. When we couple these units with repressive bonds, we can control the duration of collective oscillations for an intermediate period between collective fixedpoint behavior. The control mechanism works by monotonically increasing a single bifurcation parameter. Both the onset and arrest of oscillations are due to bifurcations. Depending on the coupling strength, the network topology and the tuning speed, our numerical simulations reveal a rich dynamics outofequilibrium with multiple inherent time scales, long transients towards the stationary states and interesting transient patterns like selforganized pacemakers. Zooming into the transition regime, we pursue the arrest of oscillations along the arms of spirals. We point out possible relations to the genetic network of segmentation clocks.
pacs:
05.45.a, 05.45.Xt, 05.65.+b1 Introduction
The motivation for our study of coupled dynamical units, each of which showing excitable or oscillatory behavior, comes from the segmentation clock. The segmentation clock stands for an oscillating multicellular genetic network that controls the sequential subdivision of the elongating vertebrate embryonic body axis into morphological somites (these are regularly sized cell clusters). The morphogenetic rhythm of the somitogenesis is based on repeated waves of oscillatory gene expression sweeping through the tissue. Usually it is the clock and wavefront mechanism [1] that is supposed to translate the oscillations in time of a clock into a periodic pattern in space when a wavefront moves posteriorly across the presomitic mesoderm (PSM) from the anterior and freezes the cellular oscillators as it passes by, recording their phases at the time of arrest. Questions of interest concern the relation between the patterns on the genetic and the cellular levels and the way of how a collective segmentation period rises from the ensemble of individual units. In [2] a delayed coupling theory was developed as a phenomenological mesoscopic description of the vertebrate segmentation in terms of phase oscillators to represent cyclic gene expression in the cells of the PSM. The oscillators are coupled to their neighbors with delay and have a moving boundary, describing the elongation of the embryo axis. They are equipped with a profile of natural frequencies that is designed to slow down the oscillations and stop them at the arrest front. In particular it is shown how the clock’s collective period depends on the delayed coupling. The very mechanism of arresting the oscillations at the arrest front is, however, not addressed, neither is the onset of oscillations described as a Hopf bifurcation from an underlying dynamical system (as, for example, in [3] or in [4]). Differently, in [5] and [6], the arrest of oscillations is supposed to result from an external signal (rather than from a suitably chosen frequency distribution or a systeminherent bifurcation).
In contrast, in our model of coupled genetic circuits, both the onset and the arrest of oscillations result from two bifurcations as one and the same bifurcation parameter is monotonically increased. Here, beyond a possible later application to the segmentation clock, we are interested in the more general question: Given a system of coupled units which individually can show excitable or limitcycle behavior, how can we control the duration of collective oscillations via bifurcations? The individual units are composed of a positive and a negative feedback loop. A first species activates itself in a positive feedback loop, and it also activates its own repressor, the second species , in a negative feedback loop. This motif is frequently found in neural and genetic networks in different realizations, whenever bistable units are coupled to negative feedback loops. Examples are signaling systems like the slime mold Dictyosthelium Discoideum [7], the embryonic division control system [8], or the MAPKcascade [9], and the circadian clock [10].
As it was shown in [11], this motif, considered now as an individual unit, shows three regimes of stationary behavior, when a single bifurcation parameter is monotonically increased: excitable behavior, where the system approaches a fixed point (if the perturbation from the fixed point exceeds a certain threshold, the system makes a long excursion in phase space before it returns to the fixed point), limitcycle behavior, and again excitable behavior. In [12] it was considered in different realizations and termed bistable frustrated unit (BFU). In the following we shall use this abbreviation as well.
When the BFUs are coupled with repressive bonds (being directed, undirected, frustrated or not frustrated), we observe the same gross collective behavior as the individual units show: a first regime now of collective fixedpoint (CFP)behavior, followed by collective oscillatory (CO)behavior, followed by a second regime of CFPbehavior. This sequence is observed for a broad range of parameters, of initial and boundary conditions, when a single, uniformly chosen bifurcation parameter is monotonically increased. Therefore the speed of its variation determines the duration of collective oscillations. The very type of collective transients and stationary states depends on the initial conditions, the boundary conditions, the network topology, and the strength of the repressive couplings, as we shall outline below. How much of the rich transient behavior can actually unfold depends on the speed of variation of the bifurcation parameter.
In view of the segmentation clock our model of coupled BFUs may match this system on an effective mesoscopic level of description, in which the interacting cells are described as coupled BFUs, arising from coupled positive and negative feedback loops on the underlying genetic level. Although first hints point in this direction, it is clearly beyond the scope of this paper to provide an exact mapping between our effective description and the segmentation clock and furthermore to identify the involved proteins and genes in the different realizations of the segmentation clock.
2 The model
The motif of a single bistable frustrated unit (BFU) is shown in figure 1. In this model, and are two different species. Although from a statistical physics point of view the identity of , is not important, we shall assume here that , are two different protein types. This is in line with the original motivation of this model as a coarsegrained description of some genetic circuits. The protein activates its own production (transcription in the biological language) and also the production of the protein , which in turn represses the production of . In this way, we have a selfactivating bistable unit which is coupled to a negative feedback loop. The simplest, idealized implementation of the BFU has been analyzed on a deterministic level, with protein concentrations as the only dynamical variables, and it is known to produce oscillations in a certain range of model parameters [13, 14]. In particular, the model studied in Ref. [14], which we adopt here, assumes that the protein production rates depend on the concentrations and of the two protein species as follows:
(1)  
(2) 
where is the ratio of the halflife of to that of . Here we shall focus on the case , that is when the protein has a much longer halflife than with a slow reaction on changes in , while has a fast response to changes in . The parameter sets the strength of repression of by . We shall assume , so that already a small concentration of will inhibit the production of . The parameter determines the basal expression level of . We set this parameter to anything larger than zero but much smaller than one, so that the system cannot be absorbed in the state and, simultaneously, the production rate of is small for . The Hill coefficients (powers of , on the r.h.s. of Eqs. (1)(2)) are chosen as in [14, 11]. The parameter is the maximal rate of production of for full activation () and no repression (). This will be our tunable parameter which we will use to control the behavior of our model, also when these units are coupled. This parameter seems to be also the easiest one to control in real, experimental systems [14]. The different fixedpoint and oscillatory regimes of such an individual unit are separated by subcritical Hopf bifurcations (for a definition see for example Ref. [15]) with corresponding hysteresis effects [11].
From now on, we do no longer distinguish between the species and their concentrations in the notation and write for the corresponding concentrations. We consider coupled BFUs with only repressing couplings (for simplicity, since the effect of activating couplings can be compensated by an appropriate choice of geometry, see [11]) according to
(3)  
Here labels the units. In some figures later we use , in which and denote the coordinates in both directions separately to facilitate the localization on the lattice. The parameters , , and are chosen uniformly over the grid, and , and with the same values as introduced for a single BFU, that is , , so without any finetuning. The parameter is also chosen independently of the lattice site, but as a bifurcation parameter it is varied between and or , depending on the repressive coupling strength, parameterized by , being weak or strong, respectively. Throughout the paper we call the coupling weak and strong, is chosen with the same value on all bonds and in opposite directions. For the variation of we distinguish three speeds: fast ( time units (t.u.)), gradual ( t.u.) and slow ( or t.u.). The adjacency matrix elements take values 1 if there is a directed coupling from to and zero otherwise. The adjacency matrix will be chosen to represent three lattice topologies as explained in the following.
The system (3) is integrated with the fourthorder RungeKutta algorithm with an integration step size of , which turned out to be sufficiently small.
Choice of topologies and boundary conditions
In general we choose free (f.b.c.) or periodic (p.b.c.) boundary conditions. For p.b.c. the lattice sites in both directions are identified modulo L, the linear extension. In order to facilitate the recognition of patterns we choose regular network topologies as indicated in figure 2: three type of rectangular lattices, usually chosen as squares of size with apart from some case studies for quasi onedimensional arrangements of up to units, and one hexagonal lattice, again of units with or . To have the option of multistable stationary states, we choose lattice and boundary conditions in a way that we have frustrated bonds for specific choices. The criteria for calling bonds frustrated in directed or undirected excitable media were presented in [11]. Accordingly the geometry of the HEXOClattice of figure 2d) induces no frustrated bonds, while all bonds are frustrated for the SQUOC, independently of the choice of boundary conditions, in figure 2a). An interesting alternative to the square lattice SQUOC is given by the SQUSC and SQUUDlattices, see figure 2 b) and c). Lattice SQUUD has no frustrated bonds for p.b.c. and all couplings chosen with the same strength in both directions of traversing a bond (so the bonds become effectively undirected), see figure 2 c). For f.b.c., we choose the bonds along the boundary to be directed, but undirected in the bulk, as displayed in figure 2 b). Apart from inducing frustration to certain bonds, the lattice SQUSC then has an overall nonvanishing chirality, and all elementary plaquettes have the same chirality when each bond in the bulk is traversed twice, but in opposite directions, as indicated in figure 2 b). In contrast, the SQUOC and the HEXOC have adjacent plaquettes with opposite chirality. These apparently minor features seem to be relevant for the patterns and their sensitivity to the boundary conditions.
In view of biological applications it is questionable whether multistability is a desired feature. Here we choose the different topologies with and without frustrated bonds and different chirality in order to illustrate the common features and the differences in the transient and stationary collective behavior of the whole system. The very choice of regular topologies may be even not unrealistic in view of cell assemblies in biological systems [16].
3 Results
3.1 Control of the duration of oscillations
Before we go into more detail on the dependence of the lattice topology, the coupling strength and the variation speed of the bifurcation parameter , we summarize the features which are in common to almost all realizations of our system. While we know from the dynamics of an individual unit that it shows three regimes when the bifurcation parameter is monotonically varied, it may come as a surprise that similar kind of regimes emerge as collective behavior, when the individual units are coupled and the very same control parameter is monotonically tuned. Note that it is again a single, now uniformly chosen control parameter for all that controls the duration of the intermediate regime of collective oscillations (CO). In the stationary state, the collective oscillations are synchronized in frequency, but differ by their patterns of synchronized phases. Now, for small and large values of we have collective fixedpoint behavior (CFP), in which all individuals approach a fixed point in concentration. This sequence of CFPCOCFP is seen for all lattice topologies, values of variation speed and couplings, see figure 3, apart from a few exceptions. Figure 3 shows snapshots of the values of for three topologies: SQUOC with p.b.c. (first row), SQUUD with f.b.c. (second row), and HEXOC with p.b.c. (third row) at time instants from the CFPregime (first and fourth column) and out of the COregime (second and third column). All simulations were run for a gradual variation of and at weak coupling, while the other parameters were kept fixed (, ). Here and throughout the paper, the color codes the value of : ranging from blue for zero, to red for a maximum value.
The collective fixed point and oscillatory behavior can be further characterized by the number of clusters, where we collect in one cluster all units for which () agree within the numerical accuracy. In this sense the uniformly blue panels represent a onecluster fixed point, panel e) a fourcluster fixed point, and the second and third columns indicate multicluster limit cycles with more or less regular patterns. While the displayed fixedpoint behavior corresponds to the stationary behavior at these values, the columns of the oscillatory regimes represent transient patterns upon approaching the stationary state. The snapshots were taken after 2500, 6000, 10000 and 15000 t.u. (first row), 100, 1500, 10000, and 15000 t.u. (second row), and 3000, 4000, 12000, and 15000 t.u. (third row). A typical time evolution of a set of units on a SQUUD lattice for a gradual variation of and weak coupling is shown in the movie [17]^{1}^{1}1Captions to the movies can be found in the following link captions.pdf .
The first exception is a HEXOC for strong coupling: instead of a CFP and COregime we find interchanging intervals of chaotic and periodic oscillations for small and intermediate values of (where the patterns of the periodic intervals are traveling waves), followed by a CFPregime for large as displayed in the movie [18]. The second exception is a SQUOClattice for weak coupling , where we find two additional regimes for small , such that the total sequence is given by CFPCOCFPCOCFP [19], characterized by a twocluster CFP, twocluster CO, onecluster CFP, multiplecluster CO, and onecluster CFP, respectively. The second and third additional regimes disappear for strong coupling.
For demonstrating the main control mechanism there was no need to await the approach of stationary states, before was changed. Therefore, we briefly summarize the different type of stationary states without entering details about the choice of parameters. In the CPFregimes we have timeindependent states, in which the concentrations approach fixedpoint values as onecluster (figure 4 b)), twocluster (not displayed), or multicluster fixed points (figure 4 a)), depending on the lattice topology and the initial conditions. In the COregime, the solutions are frequencysynchronized oscillations of individual units, arranged in two(figure 4 c))or multicluster distributions (figure 4 d)). In the latter case, the spiral patterns created in the transient phase, are preserved. For a HEXOC lattice and strong coupling we see chaotic time evolutions of concentrations rather than the usual sequence of regimes, cf. the discussion in section 3.2.2.
In general, for the lattices and fixed , the transient time to approach a fixed point is of the order of some hundred to thousand time units, while it can take up to 200000 time units to approach the synchronized limit cycles.
3.2 Transient patterns such as selforganized pacemakers
The type and duration of transient patterns depends on

the topology. We see spirals for all square lattices and no recognizable regular structure for the HEXOC, see below.

the boundary conditions. Rectangular waves are observed for f.b.c. and stripes for p.b.c. on the SQUSClattice.

the initial conditions. For all square lattices we observed multistability either only in an oscillatory regime (SQUOC) or both in a first CFP and COregime (SQUSC and SQUUD). In these regimes it therefore depends on the initial conditions to which attractor the system will evolve, and consequently which patterns will be generated. If we choose for a hexagonal lattice (HEXOC) all phases to be initially the same and start with in the oscillatory regime, the system oscillates with fully synchronized frequencies and identical phases. This limit cycle is, however, unstable: If we slightly perturb the state by kicking one phase by a difference of , the former collective uniform limit cycle is lost and the system evolves to a multicluster limit cycle. If we then start with in the CFPregime, where the fixed points are reached within only time units, and vary slowly or gradually, the system has reached the fixed point upon approaching the COregime. With these initial conditions the phases remain the same for the rest of the simulation. This is the reason why we do not see pattern formation for a HEXOClattice for a slow change of . In contrast, for fast and gradual change of and randomly chosen initial conditions, we do see irregular patterns in the COregime.
The large variety of transient patterns that is seen upon varying is due to the qualitative diversity of the initial conditions: while we have mainly chosen two type of initial conditions for fixed , random and or ordered ones, i.e. identical ones, perturbed at a single site, the initial conditions for varying are in general neither random nor highlyordered, but structured by the pattern that was reached before the value was increased. So it is the variety of transient patterns, which set the subsequent initial conditions, and the sensitivity to the initial conditions is in general pronounced, as we have multistable states as attractors with sufficiently large basins of attraction.
Apart from spirals, rectangular or plane waves, we see also target waves, emitted from units that act like dynamically generated pacemakers, this means in a selforganized way. The phenomenon is seen for a small range of speeds of increasing for the SQUOC and p.b.c., and for some values of fixed as a transient towards a twocluster limit cycle. So it is less generic than the generation of spirals, but remarkable as a transient.
In an ensemble of oscillators, pacemakers are those units, which entrain the phases of other oscillators to a synchronized oscillation in space and time such that the pacemakers become the centers of concentric target waves. Pacemakers here are observed for the SQUOClattice with p.b.c. as transient either when is fixed () and is weak, in the first COregime of the SQUOC, or when is slowly varied at weak coupling. This choice of speed within a small interval of allowed variations is essential for the observation. On the torus topology we then have two centers opposite to each other, one acting as a source of waves, the other one as their sink. The snapshots of figure 5 display the phase differences of oscillators with respect to the maximal phase value at the time instant of the snapshot. We see roughly circular spots of the same phase differences, located in the bulk of identical phases (blue background in figure 5 a)b)). The shrinking of the spot corresponds to a sink with respect to the front of constant phase differences, its extension to a source of the same front. In the snapshot of figure 5 c) we see a time instant, for which the two spots, corresponding to the source and the sink with maximal and minimal phase differences are visible at the same time, the corresponding threedimensional plot is shown in figure 5 d). The movie [20] shows an evolution of few periods on a torus. Altogether this collective arrangement of phases lasts over several thousand time units, before this transient approaches a twocluster COstate in the form of a chessboard pattern as the stationary state.
Rather remarkable about this phenomenon is its emergence in a set of units with a completely uniform choice of individual parameters () and coupling . Usually pacemakers are implemented with an ad hoc distinction via their natural frequencies (see, for example [21]) or some gradient in the natural frequencies of all oscillators [22]; in those cases the location of the pacemaker is predetermined by an adhoc implemented distinguished natural frequency. In contrast, our natural frequencies are induced by a completely uniform choice of parameters in the underlying dynamics, neither is the final location of the emerging spot distinguished in view of the p.b.c., but the very occurrence of a pacemaker breaks the symmetry between the oscillatory units, so that their heterogeneity must be dynamically generated. Our system, in spite of all couplings being repressive, is of the activatorinhibitor type, as the repression of a repressing bond acts effectively activating. Similar selforganized pacemakers have been identified in reactiondiffusion systems [23, 24, 25], there, however, as stable phenomena. In those systems the occurrence of selforganized pacemakers could be traced back to the vicinity of a supercritical pitchforkHopf bifurcation, in which a difference in frequencies is dynamically generated in this bifurcation. In our case, the conditions for the observation of these transient phenomena and the origin of their instability in terms of bifurcations should be further analyzed in future work. Although our selforganized pacemakers are transients, their duration over several thousand time units may be sufficient for providing a signal in the biological context before the system stabilizes to a stationary state.
3.2.1 A zoom into the transition regime
For the stationary states we observed two qualitatively different regimes with CFP and CObehavior. It is of particular interest, both from the viewpoint of statistical physics and of nonlinear dynamics, how the transition between the different collective behavior happens in detail. A detailed bifurcation analysis for a system of the size of units would be rather involved (so far we performed a detailed analysis for a system of only two mutually repressing units [26]). Here we want to adopt the perspective from statistical physics and analyze the emergence of individual fixedpoint behavior in the bulk of oscillatory behavior and vice versa, depending on the direction of change of the bifurcation parameter. So the first question is about a possible transient coexistence of both kind of behavior, similar to the coexistence of liquid and gaseous phases of a substrate for a given temperature and pressure.
Let us first consider the case of fixed , chosen out of a small interval in the transition regime () of the SQUUD and an even smaller interval for the SQUOC (). We then obtain transients up to time units that are characterized by coexisting, individually oscillating and “frozen” units. The behavior of individual units may change from first being oscillatory and then arrested or vice versa. Let us next consider a gradual change of over the regime . For this case we can pursue to a certain extent which oscillators get arrested first and how the arrest spreads over the lattice. Starting from a COregime (figure 6 a)), it is at the corners and the centers of the spirals (black spots in figure 6 b)), where the oscillators get arrested first, and the arrest then propagates along the boundaries and the black arms of the spirals (figure 6 c)), until a small island of oscillatory units is left in the bulk of fixedpoint behavior (red spot in figure 6 d)) that afterwards gets arrested as well. The values of the arrested oscillator phases are not those of the rotating unit at the time instant of arrest, as it is assumed in a clockwavefront mechanism, but of the fixedpoint value at the given value of . This is seen from figure 6 e) which shows the time evolution of five Avalues: the black curve corresponds to the upper left corner point which freezes first, the blue curve to the center of the spiral , which freezes next; in our selection of “freezing events” they are followed by a unit along the boundary (orange curve), where the freezing spreads, followed by a unit close to the center of the spiral (grey curve) and a unit somewhat in the bulk close to the upper right corner (red curve). The selection of the last three units is somewhat arbitrary. In principle it can be used for pursuing the chronological order of how the freezing spreads over the lattice. The time instants of the snapshots of figure 6 are 12600 (a), 13250 (b), 13500 (c) and 13710 (d). In particular panel c) illustrates the spreading along the arms of the spiral.
What distinguishes corners and the center of the spiral from other locations on the grid are the number of neighbors with a given phase value in the oscillatory regime. The less neighbors of a given unit can be in a state that would be “preferred” by this unit and stabilize its own state, the more sensitive it reacts upon a change of the parameter values, here of . Clearly it depends on the coupling sign and strength in competition with its individual dynamics, which neighboring states are favored by the unit. In systems of condensed matter it is the interface free energy between liquid and gaseous phases in comparison to the free energy in the bulk that determines the conditions for the survival of a nucleus of liquid in the bulk of gas, or vice versa. The nucleus has to exceed a critical size to survive and grow. Here it remains to quantify the condition for the “conversion” of CO into CFPregimes in terms of basins of attraction. Concretely, it is the basin of attraction of the fixedpoint state, first approached by the unit in the center of the spiral and at the corners, that should exceed a “critical size” to attract the value of A, and arresting its oscillations, so that other oscillators will follow. We shall not pursue this question here, as a further zoom into the neighborhood of first “freezing” units would be needed within the highdimensional and rough attractor landscape.
3.2.2 Influence of the coupling strength in the repressilatorlimit
Stronger repressive coupling accelerates the synchronization towards the pattern of the stationary state in the COregime. If we compare movies for the SQUOClattice (with p.b.c. and a gradual change of ) for weak and strong coupling, it is for strong coupling that the units synchronize much faster towards a state with two clusters of oscillators, that is a chessboard pattern of oscillator phases on the square lattices (cp. the movies [19] and [27]). Tuning towards the second CFPregime, the system evolves to a onecluster fixed point solution, so that the individual dynamics obviously dominates the coupling term even for strong couplings.
For strong coupling and the HEXOC topology, the lower CFPregime is completely absent and the former COregime becomes chaotic with intervals of oscillating behavior, as is best seen from the time evolution of oscillator phases, one of which is displayed in figure 7. For large enough, the system evolves to a fixed point, for the simulation displayed in figure 7 it happens at which is reached at time 80000 t.u.. For smaller values of , the system is in the limit, in which the coupling term dominates over the individual dynamics, a limit that we call repressilatorlimit.
We are particularly interested in this limit by two reasons. On the experimental side no selfactivating components have been identified yet in connection with the segmentation clock. From the theoretical perspective the individual dynamics of our BFUs is already rich and flexible, and the question arises as to whether a basic unit with a simpler individual dynamics would do it as well: when coupled to other such units, providing a simple mechanism for controlling the duration of oscillations along with a rich spectrum of transients, some of which may be relevant in the biological context as well. Therefore we studied our system in the limit of very strong coupling, where the effect of selfactivation should be suppressed and the HEXOC lattice resembles a set of repressilators [16]. Selfactivating components are completely absent in repressilators. The role of the bifurcation parameter there is played by the repressive coupling, and a COregime exists just because of the coupling. Varying now in our system the coupling (not ) from weak to strong on a HEXOC lattice for small, but fixed to realize the repressilatorlimit, we see first a CFPregime for (weak coupling), followed by chaotic behavior for intermediate coupling strengths, and again CFPbehavior for larger couplings, see figure 8. For larger lattices, the second CFPbehavior happens only for very strong coupling ( for a lattice). This means, in the limit where we can neglect the effect of selfactivation, we see chaotic or fixedpoint behavior. Our results are in agreement with the results of [16]. The other lattice topologies do not resemble the repressilatorlattice with three repressing units coupled in an elementary loop. Therefore the positive feedback components in our system of coupled BFUs seem to be an essential ingredient for the simple control mechanism to work.
3.2.3 Influence of the variation speed and boundary conditions
In view of the control mechanism we mainly focused on a gradual change of , although realized in an nonadiabatic way. The slow and fast speeds correspond to iterated quenches with short (fast speed) or longer (slow speed) time intervals of constant in between. The speed in these cases mainly influences the dynamics via the initial conditions for the subsequent value of , provided by the concentration values at the end of the preceding interval of fixed . In the case of multistability the subsequent states sensitively depend on the initial conditions. Figure 9 a) and b) illustrate the dependence on the initial conditions for otherwise the same gradual speed, where a different number of spirals is generated, for the case of a SQUUDlattice. For a fast variation of , where we change every 100 t.u., which is of the same order as the phases period, the system has not enough time to form distinguishable spirals, cf. figure 9 d) in contrast to c). In general, for a fast change of , some patterns have not enough time to evolve at all.
The influence of p.b.c. seems to be more restrictive for lattices with plaquettes of the same chirality (SQUSC figure 10 a) and c)) with p.b.c. and f.b.c., respectively, or with no chirality (SQUUD figure 10 b)) than for lattices with opposite chirality like the SQUOClattice with p.b.c. (d) and f.b.c. (e). In figure 10 a) the spirals are replaced by stripes which propagate from the boundaries and delete the spirals when coalescing with them in the bulk. Obviously the very choice of boundary conditions may have a selective effect on the evolving transient patterns that must not be ignored in view of applications to biological systems. From simulations of larger lattices of size we conclude that their influence remains pronounced for larger sizes.
3.2.4 A quasi onedimensional arrangement
In view of our original motivation in connection with the segmentation clock we also checked the type of patterns which evolve for a quasi onedimensional topology chosen as a SQUUD lattice with f.b.c.. Figure 11a) shows a multicluster fixed point, b) a multicluster oscillatory state with stripes that are traveling from the boundaries to the center of the lattice, and c) a onecluster fixed point for large . It should be noticed that the mechanism of generating the stripes here is different from the clock and wavefront mechanism of [1] that is usually supposed to generate the stripes in connection with the segmentation clock.
4 Summary and Conclusions
Coming back to our original motivation from the segmentation clock, it is two bifurcations in our system of coupled BFUs that initiate and arrest the oscillations when a single bifurcation parameter in monotonically increased. The speed of its variation determines the duration of the overall oscillation period and of the type of transient patterns. There is no need for finetuning the frequency profile via a spacedependent choice of input parameters, neither is there a need for an external field to induce an arrest of the oscillations, or for a higherlevel counting mechanism implemented in the elongating embryo. (For a careful and detailed discussion of alternative mechanisms in connection with the segmentation clock we refer to [28].) Our mechanism of pattern generation is therefore different from the clock and wavefront mechanism [1]. We can create periodic patterns in space due to the traveling waves in the oscillatory period or as a collective multicluster fixed point solution of the BFUs in the smallregime, characterized by coexisting different but frozen values of the concentrations and . Due to the repulsive coupling no neighboring oscillators share the same value unless the individual dynamics dominates the coupling term. Therefore the question arises whether there is any experimental support for our type of mechanism being at work in the segmentation clock. Given a system of coupled BFUs, each composed of a negative and positive feedback loop, the driving force for the onset and arrest of oscillations is the variation of a single bifurcation parameter. On the experimental side, retinoic acid (a morphogene) is an example of a molecule that increases in concentration monotonically when cells traverse the presomitic mesoderm from posterior to anterior, although the influence of this morphogene on oscillations is currently not (yet) understood on a molecular level and its relation to our parameter is therefore open.
The second question concerns the realization of coupled BFUs in the segmentation clock. So far, negative feedback loops were identified in the zebrafish, where certain dimers repress the genes of their components [29], but no indications for selfactivating loops have been found so far. This does, of course, not exclude their existence. On the modeling side, for comparison we simulated both, a system of coupled repressilators, and our system of coupled BFUs in the limit of very strong coupling. As outlined above, the control of the duration of oscillations failed in the sense that the first and third fixedpoint regime may be absent or the oscillatory regime is replaced by a chaotic one. So our mechanism does not seem to work with these simpler building blocks.
From a more general perspective the role of positive feedback loops in combination with negative ones was studied in [30]. There it was shown that in general a combination of interlinked positive and negative feedback loops of interacting genes allows to adjust the frequency of a negative feedback oscillator while keeping the amplitude of oscillations nearly constant. At the same time, the combination with a positive feedback loop makes the system more robust against perturbations, and therefore improves its reliability. Although it seems to be open whether these observations are relevant for the segmentation clock, these results together with our observations may trigger a search for selfactivating loops on the experimental side of the segmentation clock.
From the mere physics’ perspective our system of coupled BFUs provides a rich repertoire of transient patterns, in particular of selforganized pacemakers, and multiple inherent dynamically generated time scales with a variety of attractors. A closer zoom into the transition region, where the “conversion” between collective oscillatory and collective fixedpoint behavior happens, poses a challenge for further research.
Acknowledgments
We would like to thank F. Jülicher and D. Jörg (MaxPlanck Institute for the Physics of Complex Systems, Dresden) and L.G. Morelli (Universidad de Buenos Aires) for useful discussions about the segmentation clock. Financial support from the Deutsche Forschungsgemeinschaft (DFG, contract ME1332/191) is gratefully acknowledged.
References
References
 [1] Cooke J and Zeeman E C 1976 J. Theor. Biol. 58 455
 [2] Morelli L G, Ares S, Herrgen L, Schröter C, Jülicher F, and Oates A C 2009 HFSP Journal 3 55
 [3] Jensen M H, Sneppen K, and Tiana G 2003 FEBS Lett. 541 176
 [4] Goldbeter A, and Pourquié O 2008 J. Theor. Biol. 252 574
 [5] Goldbeter A, Gonze D, and Pourquié O 2007 Dev. Dyn. 236 1495
 [6] Santillán M and Machey M C 2008 PLoS ONE 3 e1561
 [7] Martiel J and Goldbeter A 1987 Biophys. J. 52 807
 [8] Pomerening J R, Kim S Y, and Ferrell J E Jr. 2005 Cell 122 565
 [9] Qiao L, Nachbar R B, Kevrekidis I G, and Shvartsman S Y 2007 PLoS Computational Biology 3 e184
 [10] Barkai N and Leibler S 2000 Nature 403 267
 [11] Kaluza P, MeyerOrtmanns H 2010 Chaos 20 043111
 [12] Garai A, Waclaw B, Nagel H, and MeyerOrtmanns H 2012 J. Stat. Mech. P01009
 [13] Guantes R and Poyatos J F 2006 PLoS Comp. Biol. 2 0188
 [14] Krishna S, Semsey S, and Jensen M 2009 Phys. Biol. 6 036009
 [15] Strogatz S H 1994 Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry and Engineering, Westview Press, Perseus Books Publishing LLC
 [16] Jensen M H, Krishna S, and Pigolotti S 2009 Phys. Rev. Lett. 103 118101
 [17] Supplemental material, mov1.mp4
 [18] Supplemental material, mov2.mp4
 [19] Supplemental material, mov3.mp4
 [20] Supplemental material, mov4.mp4
 [21] Radicchi F and MeyerOrtmanns H 2006 Phys. Rev. E 73 036218 17
 [22] Radicchi F and MeyerOrtmanns H 2006 Phys. Rev. E 74 026203 19
 [23] Stich M, Ipsen M, and Mikhailov A S 2001 Phys. Rev. Lett. 86 4406
 [24] Stich M, Ipsen M, and Mikhailov A S 2002 Physica D 171 19
 [25] Stich M, Mikhailov A S, and Kuramoto Y 2009 Phys. Rev. E 79 026110
 [26] Labavić D and MeyerOrtmanns H, A toggle switch of two bistable frustrated units, in preparation
 [27] Supplemental material, mov5.mp4
 [28] Oates A C, Morelli L G, and Ares S 2012 Development 139 (4) 625
 [29] Schröter C, Ares S, Morelli L G, Isakova A, Hens K, Soroldoni D, Gajewski M, Jülicher F, Maerkl S J, Deplancke B, and Oates A C 2012 PLoS Biology 10 (7) e1001364
 [30] Tsai T YC, Choi Y S, Ma W, Pomerening J R, Tang C, and Ferrell J E Jr 2008 Science 321 126