Abstract
Collective oscillations of cells in a population appear under diverse biological contexts. Here, we establish a set of common principles by categorising the response of individual cells against a timevarying signal. A positive intracellular signal relay of sufficient gain from participating cells is required to sustain the oscillations, together with phase matching. The two conditions yield quantitative predictions for the onset cell density and frequency in terms of measured singlecell and signal response functions. Through mathematical constructions, we show that cells that adapt to a constant stimulus fulfil the phase requirement by developing a leading phase in an active frequency window that enables celltosignal energy flow. Analysis of dynamical quorum sensing in several cellular systems with increasing biological complexity reaffirms the pivotal role of adaptation in powering oscillations in an otherwise dissipative celltocell communication channel. The physical conditions identified also apply to synthetic oscillatory systems.
Introduction
Homogeneous cell populations are able to exhibit a rich variety of organised behaviour, among them periodic oscillations. During mound formation of starved social amoebae, cyclic AMP waves guide migrating cells towards the high density region^{1,2,3,4,5}. Elongation of the vertebrate body axis proceeds with a segmentation clock^{6,7}. Multicellular pulsation has also been observed in nerve tissues^{8}, during dorsal closure in late stage drosophila embryogenesis^{9,10,11,12}, and more^{13}. In these examples, communication through chemical or mechanical signals is essential to activate quiescent cells. Dubbed dynamical quorum sensing (DQS) to emphasise the role of increased cell density in triggering the autoinduced oscillations, this class of behaviour lies outside the wellknown Kuramoto paradigm of oscillator synchronisation^{14,15}.
Interestingly, autoinduced oscillations have also been reported in situations without an apparent biological function. A case in point is otoacoustic emission (OAE), where a healthy human ear emits sound spontaneously in a silent environment^{16,17}. Anatomically, sound is generated by hair bundles, the sensory units of hair cells that detect sound with ultrahigh sensitivity^{18,19,20}. Another example is glycolytic oscillations of yeast cells which can be induced across different laboratory conditions^{21,22,23,24,25,26}. This type of phenotypic behaviour may not confer benefits to the organism, so their existence is puzzling.
Here, we consider a population of cells attempting to modulate temporal variations of the extracellular concentration of a protein or analyte, or a physical property of their environment, by responding to it. The response of a cell to the external property, or signal, can be mediated by an arbitrary intracellular biochemical network. By focusing on the frequencyresolved cellular response, we report a generic condition for collective oscillations to emerge, and show that it is satisfied when cells affect the signal in a way that adapts to slow environmental variations, i.e., cells respond to signal variation rather than to its absolute level. In particular, we prove the existence of an active frequency regime, where adaptive cells anticipate signal variation and attempt to amplify the signal. Sustained collective oscillations emerge when a cell population, beyond a critical density, communicates spontaneously through such a channel.
We provide a physical explanation of oscillations in terms of energydriven processes, with adaptive cells outputting energy in the active frequency regime upon stimulation. For mechanical signals, the energy output is directly observable as work on the environment. For chemical signals, chemical free energy is transferred during the release of molecules into the extracellular medium. Together with the measurable response of individual cells, quantitative predictions of the oscillation frequency and its dependence on cell density become possible.
The adaptive cellular response highlighted in this work is shown to underlie several known examples of DQS, and possibly glycolytic oscillations in yeast cell suspensions. The ubiquity of adaptation^{27,28,29,30,31,32,33,34,35,36,37} in biology may also explain the emergence of inadvertent oscillations. We discuss implications and predictions of this general mechanism at the end of the paper, in connection with previous experimental and modelling work.
Results
Necessary conditions for autoinduced oscillations
We begin by considering a scenario of mechanical oscillations, as illustrated in Fig. 1a. Later we will show that the same results hold for chemical oscillations. The cells are spatially close enough so that they could be regarded as under the same environment. Here, the extracellular signal \(s\) is taken to be the deformation of the mechanical environment, which is both sensed and modified by participating cells. The cell activity that affects the environment is denoted by a variable \(a\). Its dynamics is controlled by an unspecified intracellular regulatory network that responds to \(s\) through a mechanical sensor. To see how the intracellular activities might disrupt stasis in an equilibrium state of \(s\), we consider the following Langevin equation:
Here \(\gamma\) is the friction coefficient, \(F(s)\) the external force that tries to restore the physical environment, \(\xi\) the thermal noise, and the sum represents the total force created by \(N\) active cells in a unit volume, whose strength is set by \({\alpha }_{1}\ > \ 0\). In general, the cell activity depends on the past history of the signal. Upon a small change of \(s\), the average response of the activity of \(j\)th cell satisfies
with \(\langle \cdot \rangle\) and \({\langle \cdot \rangle }_{u}\) denoting noise average with and without an external timevarying signal, respectively. Without loss of generality, we set the stationary activity \({\langle {a}_{j}\rangle }_{u}\) to zero. The activity response function \({R}_{a}\) is a property of the intracellular molecular network, which can be computed for specific models^{38,39} or measured directly in singlecell experiments^{3,4,19,40}. In general, \({R}_{a}\) may depend on the ambient signal level \(s\) of the cell.
The shared signal \(s\) offers a means to synchronise the activities of cells. We derive here a matching condition for \(s\) and the \(a\)’s to enter a positive signal relay. Expressing Eq. (2) in Fourier form, we have \(\langle {\tilde{a}}_{j}(\omega )\rangle ={\tilde{R}}_{{a}_{j}}(\omega )\langle \tilde{s}(\omega )\rangle\). For weak disturbances, the restoring force in Eq. (1) can be approximated by a linear one, i.e., \(F(s)\simeq Ks\). Consequently, \(\langle \tilde{s}(\omega )\rangle ={\sum }_{j=1}^{N}{\alpha }_{1}{\tilde{R}}_{s}(\omega )\langle {\tilde{a}}_{j}(\omega )\rangle\), where
is the signal response function, with \(i\) the imaginary unit. For identical cells, these equations yield an oscillatory solution \(\tilde{a}({\omega }_{{\rm{o}}})\ \ne \ 0\) provided \(N{\alpha }_{1}{\tilde{R}}_{a}({\omega }_{{\rm{o}}}){\tilde{R}}_{s}({\omega }_{{\rm{o}}})=1\). To gain more insight, we express the two response functions in their amplitudes and phase shifts, i.e. \({\tilde{R}}_{a}\ \equiv \  {\tilde{R}}_{a} \exp (i{\phi }_{a})\) and \({\tilde{R}}_{s}\ \equiv \  {\tilde{R}}_{s} \exp (i{\phi }_{s})\). Then, the cell density \(N={N}_{{\rm{o}}}\) and the selected frequency \({\omega }_{{\rm{o}}}\) at the onset of collective oscillations are determined by,
These are essentially conditions of linear instability for the quiescent state expressed in terms of the singlecell and signal response functions, and constitute our first main result. For inhomogeneous cell populations, one simply replaces \({R}_{a}\) by its population average \({\overline{R}}_{a}\ \equiv \ {N}^{1}{\sum }_{j=1}^{N}{R}_{{a}_{j}}\).
Under the assumption of additive signal release from individual cells as expressed by Eq. (1), we now have a mathematical prediction for the onset density \({N}_{{\rm{o}}}\) and oscillation frequency \({\omega }_{{\rm{o}}}\). Let \({\alpha }_{2} \sim  {\tilde{R}}_{a}\) be the sensitivity of the cell activity against \(s\). We introduce a signal relay efficiency \(\overline{N}\ \equiv \ N{\alpha }_{1}{\alpha }_{2}\), which also sets the coupling strength of cellular activities through the signal. Oscillations start at the critical coupling strength \({\overline{N}}_{{\rm{o}}}={N}_{{\rm{o}}}{\alpha }_{1}{\alpha }_{2}\). Eq. (5) simply states that, at the selected frequency \({\omega }_{{\rm{o}}}\), signal amplification through the collective action of \({N}_{{\rm{o}}}\) cells compensates signal loss from dissipative forces acting on \(s\), e.g., friction for a mechanical signal or degradation/dilution for a chemical signal. The frequency \({\omega }_{{\rm{o}}}\) is chosen such that phase shifts incurred in the forward and reverse mediumcell transmissions match each other (Eq. (4)).
Although we mainly focus on the emergence of collective oscillations as cell density increases, oscillation death at high cell densities through a continuous transition, as observed in certain experimental systems, may also fulfil the selfconsistency condition (Eqs. (4–5)). An example is provided in Supplementary Note 6 (see also Supplementary Fig. 15), where, due to the nonlinear properties of the intracellular circuit, the cell activity becomes less responsive at elevated signal intensity. To make use of our procedure, the celldensity dependence of the linear response functions needs to be considered.
Celltosignal energy flow
Autoinduced collective oscillations must be driven by intracellular active processes. These active components of the system give a nonequilibrium character to the activity response^{41,42,43,44} and furthermore enable energy flow from the cell to the signal upon periodic stimulation, an interesting physical phenomenon left unnoticed so far.
To set the stage, we turn to basic considerations of nonequilibrium thermodynamics^{45,46}. The shared signal \(s\) as illustrated in Fig. 1 typically follows a dissipative dynamics such as Eq. (1). When the medium is close to thermal equilibrium, the FluctuationDissipation Theorem (FDT) relates the imaginary component \({\tilde{R}}{}_{s}^{{\prime\prime} }\) of the signal response \({\tilde{R}}_{s}\) to its spontaneous fluctuation \({\tilde{C}}_{s}\) induced by thermal noise^{38,39,47}: \(2T{\tilde{R}}{}_{s}^{{\prime\prime} }(\omega )=\omega {\tilde{C}}_{s}(\omega )\), where \({\tilde{C}}_{s}(\omega )={\langle  \tilde{s}(\omega ){ }^{2}\rangle }_{u}\) is the spectral amplitude of the signal, and \(T\) is the temperature. This relation demands \({\tilde{R}}{}_{s}^{{\prime\prime} }(\omega )\) to be positive at all frequencies. Hence, the dissipative nature of the physical environment translates into a phase delay, i.e., \({\phi }_{s}\ \equiv \ \arg ({\tilde{R}}_{s})\in (\pi ,0)\). Under the overdamped signal dynamics (Eq. (1)), Eq. (3) gives
where \({\tau }_{s}=\gamma /K\) is the signal relaxation time. (The situation \(\pi \ <\ {\phi }_{s}(\omega )\ <\ \pi /2\) occurs at high frequencies when the dynamics of \(s\) is underdamped.) On the other hand, a leading phase as required by Eq. (4) for the intracellular signal relay violates the FDT. In the present case, active cells play the role of the outofequilibrium partner. We have calculated the work done by one of the cells on the signal when the latter oscillates at a frequency \(\omega\) (see Supplementary Note 1). The output power \(\dot{W}\ \equiv \ \langle \dot{s}\cdot {\alpha }_{1}a\rangle\), i.e., the averaged value of the product between signal velocity (\(\dot{s}\)) and force from an individual cell (\({\alpha }_{1}a\)), is given by
The energy flux is positive, i.e., flowing from the cell to the signal, when \(a\) has a phase lead over \(s\), reaffirming Eq. (4) as a necessary condition on thermodynamic grounds. Stimulated energy release from an active cell to the signal as expressed by Eq. (7) constitutes our second main result in this paper.
Eq. (7) can also be used to calculate the energy flux for an arbitrary signal time series \(s(t)\), provided the linear response formula Eq. (2) applies. In particular, thermal fluctuations of \(s\) in the quiescent state may activate a net celltosignal energy flow. The total power is obtained by integrating contributions from all frequencies. Previous experiments from Hudspeth lab yielded a phaseleading response of hair bundles to mechanical stimulation at low frequencies^{19}. The same group also showed that energy can be extracted from the hair bundle via a slowly oscillating stimulus^{18}.
Chemical oscillations
The criteria given by Eqs. (4–5) apply equally to chemical oscillations illustrated in Fig. 1b. In contrast to the mechanical system, Eq. (1) at \(\gamma =1\) becomes a rate equation for the extracellular concentration \(s\) of the signalling molecules. The term \(F(s)\) (negative) gives the degradation or dilution rate of \(s\) in the medium, while individual cells secrete the molecules at a rate proportional to their activity \(a\). As the signalling molecules are constantly produced and degraded, chemical equilibrium is often violated even in the steady state. Nevertheless, \(F(s)\) usually plays the role of a stabilising force so that the signal response function \({\tilde{R}}_{s}(\omega )\) has the same phaselag behaviour as the mechanical case. Release of the molecules by the communicating cells must be phaseleading so as to drive oscillatory signalling.
Adaptive cells show phaseleading response
Apart from the aforementioned hair bundles, phaseleading response to a low frequency signal has also been reported in the activity of E. coli chemoreceptors^{40} and in the osmoresponse in yeast^{36}. Interestingly, all three of these cases are examples of adaptive sensory systems whose response to a step signal at \(t=0\) is shown in Fig. 2a. The small activity shift \(\epsilon\) at long times is known as the adaptation error. Figure 2b shows the response of the same system under a sinusoidal signal. The low frequency response exhibits a phase lead while the high frequency one has a phase lag. Below, we show that the sign switch in the phase shift of an adaptive variable is an inevitable consequence of causality.
From the causality condition \({R}_{a}(t\ <\ 0)=0\), the real (\({\tilde{R}}{}_{a}^{{\prime} }\)) and imaginary (\({\tilde{R}}{}_{a}^{{\prime\prime} }\)) part of the response function in frequency space satisfy the KramersKrönig relation^{48}:
For a step signal of unit strength, Eq. (2) yields
Comparing Eqs. (8) and (9) in the limit \(\omega \to 0\) and assuming \(\epsilon\) to be sufficiently small, we see that \({\tilde{R}}{}_{a}^{{\prime\prime} }(\omega )\) inside the integral must change sign. In other words, both phaseleading (\({\tilde{R}}{}_{a}^{{\prime\prime} }\ <\ 0\)) and lagged (\({\tilde{R}}{}_{a}^{{\prime\prime} }\ > \ 0\)) behaviour are present across the frequency domain. This is our third result.
Adaptation plays a key role in biochemical networks^{27,28}, and especially in sensory systems^{20,29,31,32,33,34,35,36}. Connection between adaptation and collective oscillations has been implicated in previous works^{4,49,50}. With the mathematical results presented above, the pathway from adaptation to phaseleading response, and onto collective oscillations through signal relay, is firmly established (See “Methods”). Below we illustrate, with the help of three examples of increasing complexity, how this line of reasoning could link up different aspects of system behaviour to arrive at a renewed understanding. Implications of our model study to experimental work are given in the “Discussion” section.
A weakly nonlinear model with adaptation
We consider first a noisy twocomponent circuit which is a variant of the model for sensory adaptation in E. coli^{41} (Fig. 3a, see also “Methods”). For weak noise, the intracellular signal relay in the quiescent state is essentially linear with the receptor response function given by
where
Here \({\tau }_{a}\) and \({\tau }_{y}\) are the timescales for the activity (\(a\)) and negative feedback (\(y\)) dynamics, respectively. In Figs. 3b, c, we show the phase shift \({\phi }_{a}(\omega )\) and the real and imaginary part of \({\tilde{R}}_{a}(\omega )\) against the frequency \(\omega\), plotted on semilog scale. As predicted, \({\phi }_{a}(\omega )\) undergoes a sign change at \({\omega }^{* }\). Correspondingly, the imaginary component of the response \({\tilde{R}}{}_{a}^{{\prime\prime} }\) becomes negative in the phaseleading regime, violating the FDT. The peak of \( {\tilde{R}}_{a}(\omega )\) is located close to \({\omega }^{* }\), with a relative width \(\Delta \omega /{\omega }^{* }\simeq {Q}^{1}\) where \(Q={\tau }_{a}{\omega }^{* }\simeq {({\tau }_{a}/{\tau }_{y})}^{1/2}\).
Allowing the chemoreceptor activity \(a\) to affect the signal as in Eq. (1) with \(F(s)=Ks\), we observe an oscillatory phase upon increase in cell density in numerical simulations (Fig. 3d). Figure 3e shows the oscillation amplitude (upper panel) and frequency (lower panel) against the coupling strength \(\overline{N}\) around the onset of oscillations. The threshold coupling strength \({\overline{N}}_{{\rm{o}}}={N}_{{\rm{o}}}{\alpha }_{1}{\alpha }_{2}\) and the onset frequency \({\omega }_{{\rm{o}}}\) both agree well with the values predicted by Eqs. (4)–(5) (see arrows in Fig. 3e). The transition is well described by a supercritical Hopf bifurcation. At finite oscillation amplitudes, there is a downward shift of the oscillation frequency which can be quantitatively calculated in the present case by introducing a renormalised response function \({\tilde{R}}_{a}^{+}(\omega )\) (see Supplementary Note 2), whose phase is shown in Fig. 3f. The oscillation frequency is determined by the crossing of the two curves \({\phi }_{s}(\omega )\) and \({\phi }_{a}^{+}(\omega )\), with the formal independent of the oscillation amplitude \(A\). As the oscillation amplitude grows further, higher order harmonics generated by the nonlinear term become more prominent. The model system eventually exits from the limit cycle through an infiniteperiod bifurcation and arrives at a new quiescent state. The upper bifurcation point \({\overline{N}}_{b}\) is inversely proportional to the adaptation error \(\epsilon\) (see Supplementary Fig. 1).
The signal phase shift \({\phi }_{s}(\omega )\) is given by Eq. (6). When the signal relaxation time \({\tau }_{s}\) is much shorter than the cell adaptation time \({\tau }^{* }\ \equiv \ 2\pi /{\omega }^{* }\), \({\phi }_{s}(\omega )\) stays close to zero so that the selected period is essentially given by \({\tau }^{* }\). In this case, \( {\tilde{R}}_{a}(\omega )\) is near its peak and hence the cell density required by Eq. (5) is the lowest. As signal clearance slows down, the crossing point shifts to lower frequencies. Given a finite adaptation error \(\epsilon \ > \ 0\), there is a generic maximum signal relaxation time \({\tau }_{s}^{* } \sim {\epsilon }^{1}\) beyond which the phase matching cannot be achieved (see Supplementary Note 3 and also Supplementary Fig. 3).
Excitable dynamics
DQS in Dictyostelium and other eukaryotic cells takes the form of pulsed release of signalling molecules^{2,7,51}. The highly nonlinear twocomponent FitzHughNagumo (FHN) model is often employed for such excitable phenomena^{3,52,53}. Similar to the sensory adaptation model discussed above, each FHN circuit has a memory node \(y\) that keeps its activity \(a\) low (the resting state) under a slowvarying signal \(s(t)\) (Fig. 4a, see also “Methods”). On the other hand, a sufficiently strong noise fluctuation or a sudden shift of \(s\) sends the circuit through a large excursion in phase space (known as a firing event) when \(y\) is slow (i.e., \({\tau }_{y}\gg {\tau }_{a}\)). Our numerical investigations show that the noisetriggered firing does not disrupt the adaptive nature of the circuit under the negative feedback from \(y\). The noiseaveraged response of a single FHN circuit exhibits the same characteristics as the sensory adaptation model, including adaptation to a stepwise stimulus after a transient response (Fig. 4b, upper panel), as well as the phaseleading behaviour and diminishing response amplitude on the low frequency side (Fig. 4b, lower panel).
Fig. 4c shows time traces of individual cell activities (blue and green curves) as well as that of the signal \(s\) (red curve) from simulations of weakly coupled FHN circuits at three different values of the coupling strength \(\overline{N}\) (see Methods). At \(\overline{N}=0.5\), the two selected cells fire asynchronously while \(s\) remains constant. At \(\overline{N}=0.9\), collective behaviour as seen in the oscillation of \(s\) starts to emerge, although individual circuits continue to fire sporadically. Upon further increase of \(\overline{N}\), synchronised firing is established. Despite the highly nonlinear nature of the FHN model, both the onset coupling strength \({\overline{N}}_{{\rm{o}}}\) and the frequency \({\omega }_{{\rm{o}}}\) are well predicted by Eqs. (4)–(5) using the respective response functions in the resting state (Fig. 4d).
Yeast glycolytic oscillations
We take the adapttooscillate scenario one step further to examine the dynamics of ATP autocatalysis in yeast. Concentration oscillations of NADH and glycolytic intermediates have been observed in yeast cell extracts as well as in starved yeast cell suspensions upon shutting down the respiratory pathway (see ref. ^{22} for a review). The phosphofructokinase (PFK), an enzyme in the upper part of the glycolytic pathway, is tightly regulated by ATP, a key product of glycolysis. This robust negative feedback is commonly regarded as the driver of glycolytic oscillations, with a typical period of 30–40 s in intact cells but 2 min or longer in extracts. Cells at high density oscillate synchronously due to redox signalling via the freely diffusing molecule acetaldehyde (ACE)^{22,26}. As the cell density decreases, the synchronised behaviour breaks down. While many studies found continued oscillation of individual cells at their own frequencies^{25,54}, simultaneous disappearance of individual and collective oscillations as in other DQS systems has also been reported^{23}. We show below that both types of behaviour could be accommodated in a model of glycolysis that couples intracellular redox state to ATP autocatalysis.
We first investigate the dynamic properties of singlecell glycolysis at fixed intracellular glucose and ACE concentrations. Collective oscillations in yeast cell suspensions, which require ACE transport across the cell membrane, will be discussed later. Our starting point is the du Preez et al. model^{55} that includes around 20 metabolic reactions (Fig. 5a). By monitoring the temporal response of metabolites under perturbations of the intracellular ACE concentration, we obtained a phase diagram shown in Fig. 5b. The white region marks spontaneous oscillations in an isolated cell^{56}. The glucose concentration, which controls glycolytic flux, needs to be sufficiently high for oscillations to take place. ACE also has a role in the dynamics: either very low or very high concentrations arrest the oscillations.
The nonoscillatory part of the phase diagram can be further divided into three subregimes according to the adaptive properties of the metabolic network against an ACE signal (coloured bands in Fig. 5b). Figure 5c gives, under one representative condition in each band, the concentration variation of four metabolites upon a sudden shift in the intracellular ACE concentration ACE\({}_{{\rm{in}}}\). In all cases, the intracellular redox agent NAD follows closely ACE concentration change and hence acts as an instantaneous transducer of the signal. ATP adapts in both the orange and blue band, while PYR, the substrate to produce ACE, adapts only in the blue band. TRIO, the metabolite immediately upstream of the reaction GAPDH that uses NAD as cofactor, does not adapt. Overall, the adaptation error increases progressively as one moves away from the oscillatory region. A sign change in the response of ATP (and also of TRIO) takes place across the dashed line at ACE\({}_{{\rm{in}},0}\simeq 0.2\) mM.
Figure 5d shows phase shifts of ATP, NAD and five other metabolites along the glycolytic pathway against a periodic ACE signal at various frequencies. At the point marked by star in the orange band, ATP, BPG and PEP are phaseleading (after a \(\pi\) shift) below the frequency \({\omega }^{* }\simeq 20\) min\({}^{1}\) (upper panel). The list is expanded to all six metabolites (except NAD which is synchronised with the signal) when the environment shifts to the point marked by diamond in the blue band (lower panel). In particular, PYR (blue line) have a large phase lead around \({\omega }^{* }\), which matches its adaptive behaviour under Fig. 5c (Supplementary Note 4).
To further disentangle dynamical properties of the network, we constructed a reduced model in Fig. 5e by taking into account stoichiometry and known regulatory interactions along the glycolytic pathway^{24}, and by making use of the timescale separation in the turnover of metabolites as suggested by their response spectra (Fig. 5d and Supplementary Note 5). Since ATP and PYR now appear as coproducts of the condensed reaction PYK in the reduced model, the latter can be viewed as a reporter of ATP homeostasis implemented by the negative feedback loop (cyan line in Fig. 5e). Figure 5f shows the phase diagram of the reduced model against the intracellular ACE\({}_{{\rm{in}},0}\). Similar to the full model at high glucose concentrations, the circuit enters an oscillatory state at intermediate values of ACE\({}_{{\rm{in}},0}\), and shows adaptive response on the two wings. The extended adaptive regime on the high ACE\({}_{{\rm{in}},0}\) side differs from the behaviour seen in Fig. 5b, but is reproduced by a mutant of the full model where the glyoxylate shunt (GLYO) is turned off (Supplementary Figs. 10–11).
We now examine a model of yeast cell suspensions where individual cells metabolise according to the reduced model and communicate their redox state through ACE (Supplementary Note 6). ACE is synthesised internally and degraded at rates \({k}_{{\rm{in}}}\) and \({k}_{{\rm{ex}}}\) within and outside the cell, respectively. The rate of its crossmembrane transport is set by the membrane permeability \(D\). When \(D\) is large, the intracellular and extracellular ACE levels are essentially the same (left panel, Fig. 5g). On the low density side, the homogeneous cell population enters the oscillatory phase through synchronisation of oscillatory cells. This behaviour continues beyond the point ACE\({}_{{\rm{in}}}=0.72\) (arbitrary unit) when an isolated cell switches from oscillation to adaptation (Fig. 5f). In other words, the system crosses over smoothly from oscillator synchronisation to adaptationdriven oscillations, or DQS. As \(D\) decreases and becomes comparable to ACE degradation rates (set at \({k}_{{\rm{in}}}=0.5\) and \({k}_{{\rm{ex}}}=0.3\)), the intracellular ACE concentration grows and eventually exceeds the upper selfoscillatory threshold ACE\({}_{{\rm{in}}}=0.72\) even at the low cell density limit. Our simulations indicate that DQS persists at \(D=0.4\) but disappears at \(D=0.1\) (middle and right panels, Fig. 5g). Interestingly, DQS at \(D=0.4\) disappears at an upper threshold density \({\rho }_{{\rm{c}}2}=0.73\). This inverse DQS, i.e., oscillation quenching at high cell density, can be quantitatively explained by Eqs. (4)–(5) using the numerically determined response functions that depend on the cell density (Supplementary Note 6 and Supplementary Fig. 15).
Discussion
In this work, we investigated a general scenario for emerging oscillations in a group of cells that communicate via a shared signal. It covers a broad class of pulsation behaviour in cell populations, collectively known as DQS. Using the singlecell response to external stimulation, we formulated a quantitative requirement for the onset of collective oscillations that must be satisfied by active cells as well as models of them. A proof is presented to link this requirement to the adaptive release of signalling molecules by individual cells. Our work thus consolidates observations made in the literature and formalises adaptation as a unifying theme behind DQS.
The above mathematical results connect well to the recent surge of interest in active systems, where collective phenomena emerge due to energydriven processes on the microscopic scale^{57,58}. The study of such nonequilibrium processes opens a new avenue to explore mechanisms of spontaneous motion on large scales. We presented a general formula for the energy outflow of a living cell through a designated mechanical or chemical channel under periodic stimulation. This energy flux is positive over a range of frequencies when the cell responds to the stimulus adaptively. Since adaptation is a measurable property of a cell, the thermodynamic relation is applicable without making specific assumptions about intracellular biochemical and regulatory processes, while most models do. When cells are placed together in a fixed volume, a quorum is required to activate the energy flow via self and mutual stimulation.
We reported three case studies to illustrate how these general yet quantitative relations could be applied to analyse the onset of collective oscillations in specific cell populations. Our first example is a coarsegrained model where signal reception and release are integrated into the same activity node (e.g., a membrane protein or a molecular motor). Due to the weak nonlinearity of the intracellular circuit, many analytical results were obtained. The intracellular adaptive circuit has two timescales: the activity relaxation time \({\tau }_{a}\) and the negative feedback time \({\tau }_{y}\). Their ratio \({Q}^{2}={\tau }_{a}/{\tau }_{y}\), similar to the quality factor in resonators, determines the shape of the adaptive response (Fig. 2). At small adaptation error \(\epsilon \ll 1\), the imaginary part of the response function \({\tilde{R}}_{a}(\omega )\) changes sign at the characteristic frequency \({\omega }^{* }\simeq {({\tau }_{a}{\tau }_{y})}^{1/2}\). This is also approximately the frequency where \( {\tilde{R}}_{a}(\omega )\) reaches its maximum. When cells are coupled through the signal with a relaxation time \({\tau }_{s}\), the onset oscillation frequency \({\omega }_{{\rm{o}}}\) increases with decreasing \({\tau }_{s}\), reaching its maximal value \({\omega }^{* }\) when \({\tau }_{s}\ll 1/{\omega }^{* }\).
Much of these results carry over to our second example, a population of coupled excitable circuits described by the FHN model. Despite its highly nonlinear nature, the FHN model in the resting state shows adaptive response under weak stimulation. Our numerical simulations of the coupled system at weak noise confirm the onset oscillation frequency and the critical cell density predicted by Eqs. (4)–(5).
The above theoretical predictions compare favourably with available experimental data. The first is mechanical stimulation of hair cells carried out by Martin et al.^{19}, where the cellular response was extracted using a flexible glass fibre. Deformation of the glass fibre, which is the signal here, has a relaxation timescale (\(\sim\! 0.5\) ms) much shorter than the adaptation time of the hair bundle (\(\sim\! 0.1\) s). Spontaneous oscillations of the combined system were observed at 8 Hz, the predicted frequency where the imaginary part of the hair bundle response function \({\tilde{R}}{}_{a}^{{\prime\prime} }(\omega )\) undergoes the expected sign change. The second is a recent microfluidic singlecell measurement of Dictyostelium reported by Sgro et al.^{3}, where the change of cytosolic cAMP level (activity \(a\)) in response to extracellular cAMP variation (signal \(s\)) was presented. From the measured response \(a(t)\) to a step increase of the signal in their work (reproduced in Fig. 6a, upper panel), we computationally deduced the response function \({R}_{a}(t)={\mathrm{d}}a/\mathrm{d}t\) in the time domain (Fig. 6a, lower panel) and then the response spectrum \({\tilde{R}}_{a}\) via Fourier transform. The resulting phase shift \({\phi }_{a}\) changes sign around \({\omega }^{* }=1\) min\({}^{1}\) (Fig. 6b). According to our theory, the onset oscillation period at high flow rate should be around 6.28 min, which is indeed what was observed in experiments^{2,3,4}.
DQS in Dictyostelium is a timedependent phenomenon coupled to cell migration and development^{1,5}. In the experiments reported in refs. ^{2,3}, synchronised firing of cells starts five hours after nutrient deprivation. The period of firing shortens from 15–30 min at the onset to 8 min and thereafter 6 min as cells begin to aggregate. Therefore the onset of collective oscillations may not be triggered by a critical cell density as such but the cell density does affect the period of oscillations. Previously, a property known as foldchange detection (FCD) was invoked and verified to explain cellcell signalling even when cells are far apart^{4}. In FCD, the intracellular signal relay circuit is activated by a relative change \(\Delta s/s\) of the signal \(s\). Consequently, the detection sensitivity of the activity response function \({\alpha }_{2} \sim 1/s\). In a population of communicating cells, the signal strength \(s\) is proportional to the cell density \(N\). Hence, FCD renders the signal relay efficiency \(\overline{N}\) independent of the cell density \(N\). To explain the accelerated pulsing at increasing cell densities, other aspects of the system need to be considered, e.g., cAMP clearance by phosphodiesterase secreted by cells^{61}. Building these details into the FHN model, Sgro et al.^{3} showed that the coupled equations are able to qualitatively reproduce the observed behaviour. The data analysis procedure illustrated by Fig. 6 offers a direct way to link pulsation from 8 min to about 6 min with a faster signal clearance effected by a higher concentration of phosphodiesterase in the surrounding medium. The long firing interval at early stage of the development could be attributed to physiological differences in the intracellular molecular network, e.g., a much longer negative feedback time \({\tau }_{y}\) that awaits experimental verification^{59,60}. With this type of data, a similar procedure could be applied to analyse the segmentation clock in the presomitic mesoderm^{7}.
Our third example, the glycolytic oscillation in yeast cell suspensions, is also an open problem. Simulation studies of a detailed model of yeast glycolysis^{55} yielded a relatively simple phase diagram shown in Fig. 5b, with the intracellular glucose and ACE concentrations as control parameters. As reported previously^{55}, cells in the white region oscillate spontaneously in a constant environment, driven by an instability associated with the negative feedback in ATP autocatalysis. In the neighbourhood of this region, we found that the ATP concentration adapts to the intracellular environment, in particular to a sudden shift in ACE concentration that affects directly the intracellular NAD/NADH ratio. The adaptation error increases as one moves away from the oscillatory region. These dynamical features are captured by a reduced model of ATP autocatalysis we proposed to approximate the lowdimensional attractor of the full model at high glucose concentrations. We then considered a homogeneous population of cells that carry out fermentation according to the reduced model, using membrane permeability \(D\) to tune intracellular ACE concentration at a given cell density. When \(D\) is much greater than the ACE turnover/degradation rates, the intracellular and extracellular ACE levels are equilibrated. In such a situation, collective oscillations on the low cell density side first emerge through synchronisation of individual cells that enter the selfoscillatory state. Further increase of the cell density elevates both intracellular and extracellular ACE levels, eventually brings individual cells out of the selfoscillatory state. However, the population continues to oscillate following the DQS scenario. Slower crossmembrane diffusion drives up intracellular ACE level and, at some point, eliminates selfoscillation. Nevertheless, the population may still oscillate via DQS at an intermediate range of cell densities. Beyond an upper critical cell density, the diminishing adaptive response of glycolytic flux to the signal eventually arrests collective oscillations. This new phenomenon, which we name inverse DQS, is quantitatively predicted by our theory. We note that the strong celltocell variability observed in singlecell experiments^{56} could significant alter the behaviour shown in Fig. 5g on the low cell density side, an issue we leave to future work.
These model studies helped to refine and resolve various quantitative issues in the induction of collective oscillations in wellstudied systems, and at the same time inspire novel applications built around adaptationdriven signal relay. One promising direction to follow is the development of artificial oscillatory systems with techniques from synthetic biology^{62,63,64,65}. In analogy with the hair cell/glass fibre setup, one may think of tricking a quorumsensing cell to oscillate by confining it to a volume small enough to enable positive signal relay.
In statistical physics, the response function formalism is widely used to analyse system level response to environmental perturbations, but its application to collective behaviour in biological systems is still limited. Our examples show that cell models with different levels of biological detail, out of either necessity or convenience, could yield qualitatively or even quantitatively similar response curves with respect to, say the production of a particular chemical used in celltocell communication, which is reassuring. As these curves are increasingly accessible from experiments, their direct use for analysis and hypothesis building is highly desirable. With respect to the link between adaptation and collective oscillations, our formulation unifies and generalises previous studies in at least three specific settings. The first is an abstract 3variable model that connects foldchange detection of individual cells to the robustness of collective oscillations over a broad range of cell densities^{4}. In the second case, adaptation was proposed to play an important role in the collective oscillation of neuronal networks^{49}. Lastly, an Isingtype model of chemoreceptor arrays in E. coli^{50} predicts that increasing the coupling strength between adaptive receptors drives the system to collective oscillations, although in reality the chemoreceptor array manages to operate below the oscillatory regime. Despite the risk of running into an oscillatory instability, the coupling enhances sensitivity of the array to ligand binding. Along this sensitivitystability tradeoff, one may speculate that some of the reported collective oscillations under laboratory conditions could actually arise from over perfection of adaptive/homeostatic response in the natural environment, a hypothesis that invites further experimental testing.
Methods
An adaptive model with cubic nonlinearity
The data presented in Fig. 3 were obtained from numerical integration of the coupled equations^{41,44}: \({\tau }_{a}\dot{a}=a{c}_{3}{a}^{3}+y+{\alpha }_{2}s+{\eta }_{a},\) and \({\tau }_{y}\dot{y}=a\epsilon y+{\eta }_{y}.\) Here \(y\) is a memory node that implements negative feedback control on \(a\), \(\epsilon\) sets the adaptation error, and \({\tau }_{a}\) and \({\tau }_{y}\) are the intrinsic timescales for the dynamics of \(a\) and \(y\), respectively. \({\eta }_{a}\) and \({\eta }_{y}\) are gaussian white noise with zero mean and correlations: \(\langle {\eta }_{a}(t){\eta }_{a}(\tau )\rangle =2T{\tau }_{a}\delta (t\tau )\) and \(\langle {\eta }_{y}(t){\eta }_{y}(\tau )\rangle =2T{\tau }_{y}\delta (t\tau )\), where \(\delta (t)\) is the Dirac delta function. The cubic nonlinearity (\({c}_{3}{a}^{3}\)) is needed to limit cellular activity to a finite strength. For simplicity, we choose \({\alpha }_{2}=1\) so that the response function defined by \({\tilde{R}}_{a}(\omega )=\langle \tilde{a}(\omega )\rangle /\tilde{s}(\omega )\) can be compared with its equilibrium counterpart that satisfies the FDT \({\tilde{R}}{}_{a}^{{\prime\prime} }=\omega {\tilde{C}}_{a}(\omega )/(2T)\), with \({\tilde{R}}{}_{a}^{{\prime\prime} }\) denoting the imaginary component of \({\tilde{R}}_{a}\). Data in Fig. 3 were obtained by coupling cells via Eq. (1) with \(F(s)=Ks\) and \(\xi =0\).
Existence of oscillatory state under an adaptive response
We have shown in the Main Text that adaptive intracellular observables exhibit a phaseleading response in a certain frequency interval. For a given adaptive observable \(a\), the phase lead \({\phi }_{a}(\omega )\) spans a continuous range from 0 to a maximum value \({\phi }_{a}^{\max }\) (\(<\ \pi\)). Meanwhile, the phase delay \({\phi }_{s}={\tan }^{1}(\omega {\tau }_{s})\) varies continuously from 0 to \(\pi /2\) (Eq. (6)). Since \({\tau }_{s}\) controls how fast \({\phi }_{s}(\omega )\) decreases from 0 to \(\pi /2\) as \(\omega\) increases, intersection of \({\phi }_{s}(\omega )\) with \({\phi }_{a}(\omega )\) can always be found by tuning \({\tau }_{s}\). In particular, when \({\tau }_{s}\to 0\), a solution is found at the high frequency end of the active frequency interval where \({\phi }_{a}(\omega )=0\). From this discussion, we see that the onset frequency \({\omega }_{{\rm{o}}}\) of oscillations is mostly determined by the intracellular dynamics, i.e., \({\phi }_{a}(\omega )\), but the medium can have a weak effect on \({\omega }_{{\rm{o}}}\) when its relaxation time is comparable to that of the intracellular dynamics.
Coupled FHN model
A single FHN circuit takes the form, \({\tau }_{a}{\dot{a}}_{j}={a}_{j}{a}_{j}^{3}/3{y}_{j}+{\alpha }_{2}s+{\eta }_{{a}_{j}}\), \({\tau }_{y}{\dot{y}}_{j}={a}_{j}\epsilon {y}_{j}+{a}_{0}+{\eta }_{{y}_{j}}.\) The positive sign of the first term in the equation for \({a}_{j}\) gives rise to excitability. In the absence of the stimulus \(s\), each cell assumes the resting state with a mean activity \({a}_{rs}\ \equiv \ \langle {a}_{j}(t)\rangle\). For small values of \(\epsilon\), the resting state activity \({a}_{rs}\simeq {a}_{0}\) is nearly constant under a slowvarying \(s(t)\). FHN circuits are coupled together through a signal field whose dynamics is described by, \({\tau }_{s}\dot{s}=s+{\alpha }_{1}{\sum }_{j}^{N}({a}_{j}{a}_{rs})\). The parameters used in generating Fig. 4 are: \({\alpha }_{2}=1\), \(N=1000\), \(\epsilon =0.1\), \(T=0.1\), \({\tau }_{a}=1\), \({\tau }_{y}=5\), \({a}_{0}=1.5\), and \({\tau }_{s}=1\). \({\alpha }_{1}=\overline{N}/(N{\alpha }_{2})\) is determined by the control parameter \(\overline{N}\).
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Data availability
The data that support the findings of this study are available from S.W. Wang on request. They can also be generated from the provided code.
Code availability
The code that support the findings of this study are available at https://github.com/ascendancy09/Collectiveoscillations.
References
 1.
Schaap, P. Evolutionary crossroads in developmental biology: dictyostelium discoideum. Development 138, 387–396 (2011).
 2.
Gregor, T., Fujimoto, K., Masaki, N. & Sawai, S. The onset of collective behavior in social amoebae. Science 328, 1021–1025 (2010).
 3.
Sgro, A. E. et al. From intracellular signaling to population oscillations: bridging sizeand timescales in collective behavior. Mol. Syst. Biol. 11, 779 (2015).
 4.
Kamino, K. et al. Foldchange detection and scale invariance of cell–cell signaling in social amoeba. Proc. Natl Acad. Sci. USA 114, 201702181 (2017).
 5.
Bretschneider, T., Othmer, H. G. & Weijer, C. J. Progress and perspectives in signal transduction, actin dynamics, and movement at the cell and tissue level: lessons from dictyostelium. Interface Focus 6, 20160047 (2016).
 6.
Hubaud, A. & Pourquié, O. Signalling dynamics in vertebrate segmentation. Nat. Rev. Mol. Cell Biol. 15, 709 (2014).
 7.
Hubaud, A., Regev, I., Mahadevan, L. & Pourquie, O. Excitable dynamics and yapdependent mechanical cues drive the segmentation clock. Cell 171, 668–682 (2017).
 8.
Buzsaki, G. Rhythms of the Brain (Oxford University Press, 2006).
 9.
Solon, J., KayaÇopur, A., Colombelli, J. & Brunner, D. Pulsed forces timed by a ratchetlike mechanism drive directed tissue movement during dorsal closure. Cell 137, 1331–1342 (2009).
 10.
Sokolow, A., Toyama, Y., Kiehart, D. P. & Edwards, G. S. Cell ingression and apical shape oscillations during dorsal closure in drosophila. Biophys. J. 102, 969–979 (2012).
 11.
Dierkes, K., Sumi, A., Solon, J. & Salbreux, G. Spontaneous oscillations of elastic contractile materials with turnover. Phys. Rev. Lett. 113, 148102 (2014).
 12.
Jülicher, F. & Prost, J. Spontaneous oscillations of collective molecular motors. Phys. Rev. Lett. 78, 4510 (1997).
 13.
Ko, C. H. et al. Emergence of noiseinduced oscillations in the central circadian pacemaker. PLoS Biol. 8, e1000513 (2010).
 14.
Kuramoto, Y. Chemical Oscillations, Waves, and Turbulence, vol. 19 (Springer Science & Business Media, 2012).
 15.
Strogatz, S. Sync: The Emerging Science of Spontaneous Order (Penguin UK, 2004).
 16.
Gold, T. Hearing. ii. the physical basis of the action of the cochlea. Proc. R. Soc. Lond. B 135, 492–498 (1948).
 17.
Kemp, D. T. Stimulated acoustic emissions from within the human auditory system. J. Acoust. Soc. Am. 64, 1386–1391 (1978).
 18.
Martin, P. & Hudspeth, A. J. Active hairbundle movements can amplify a hair cell’s response to oscillatory mechanical stimuli. Proc. Natl Acad. Sci. USA 96, 14306–14311 (1999).
 19.
Martin, P., Hudspeth, A. & Jülicher, F. Comparison of a hair bundle’s spontaneous oscillations with its response to mechanical stimulation reveals the underlying active process. Proc. Natl Acad. Sci. USA 98, 14380–14385 (2001).
 20.
Hudspeth, A. Integrating the active process of hair cells with cochlear function. Nat. Rev. Neurosci. 15, 600 (2014).
 21.
Goldbeter, A. Biochemical Oscillations and Cellular Rhythms (Cambridge University Press, 1997).
 22.
Richard, P. The rhythm of yeast. FEMS Microbiol. Rev. 27, 547–557 (2003).
 23.
De Monte, S., d’Ovidio, F., Danø, S. & Sørensen, P. G. Dynamical quorum sensing: Population density encoded in cellular dynamics. Proc. Natl Acad. Sci. USA 104, 18377–18381 (2007).
 24.
Chandra, F. A., Buzi, G. & Doyle, J. C. Glycolytic oscillations and limits on robust efficiency. Science 333, 187–192 (2011).
 25.
Gustavsson, A.K. et al. Sustained glycolytic oscillations in individual isolated yeast cells. FEBS J. 279, 2837–2847 (2012).
 26.
Amemiya, T. et al. Collective and individual glycolytic oscillations in yeast cells encapsulated in alginate microparticles. Chaos 25, 064606 (2015).
 27.
Goentoro, L. & Kirschner, M. W. Evidence that foldchange, and not absolute level, of \(\beta\)catenin dictates wnt signaling. Mol. Cell 36, 872–884 (2009).
 28.
CohenSaidon, C., Cohen, A. A., Sigal, A., Liron, Y. & Alon, U. Dynamics and variability of erk2 response to egf in individual living cells. Mol. Cell 36, 885–893 (2009).
 29.
Alon, U., Surette, M. G., Barkai, N. & Leibler, S. Robustness in bacterial chemotaxis. Nature 397, 168–171 (1999).
 30.
Tu, Y. Quantitative modeling of bacterial chemotaxis: Signal amplification and accurate adaptation. Annu. Rev. Biophys. 42, 337 (2013).
 31.
Reisert, J. & Matthews, H. R. Response properties of isolated mouse olfactory receptor cells. J. Physiol. 530, 113–122 (2001).
 32.
Hohmann, S. Osmotic stress signaling and osmoadaptation in yeasts. Microbiol. Mol. Biol. Rev. 66, 300–372 (2002).
 33.
Hazelbauer, G. L., Falke, J. J. & Parkinson, J. S. Bacterial chemoreceptors: highperformance signaling in networked arrays. Trends Biochem. Sci. 33, 9–19 (2008).
 34.
Menini, A. Calcium signalling and regulation in olfactory neurons. Curr. Opin. Neurol. 9, 419–426 (1999).
 35.
Nakatani, K., Tamura, T. & Yau, K. Light adaptation in retinal rods of the rabbit and two other nonprimate mammals. J. Gen. Physiol. 97, 413–435 (1991).
 36.
Mettetal, J. T., Muzzey, D., GómezUribe, C. & van Oudenaarden, A. The frequency dependence of osmoadaptation in Saccharomyces cerevisiae. Science 319, 482–484 (2008).
 37.
Hoeller, O., Gong, D. & Weiner, O. D. How to understand and outwit adaptation. Dev. cell 28, 607–616 (2014).
 38.
Wang, S.W., Kawaguchi, K., Sasa, S.i & Tang, L.H. Entropy production of nanosystems with time scale separation. Phys. Rev. Lett. 117, 070601 (2016).
 39.
Wang, S.W. Inferring energy dissipation from violation of the fluctuationdissipation theorem. Phys. Rev. E 97, 052125 (2018).
 40.
Shimizu, T. S., Tu, Y. & Berg, H. C. A modular gradientsensing network for chemotaxis in Escherichia coli revealed by responses to timevarying stimuli. Mol. Syst. Biol. 6, 382 (2010).
 41.
Lan, G., Sartori, P., Neumann, S., Sourjik, V. & Tu, Y. The energyspeedaccuracy tradeoff in sensory adaptation. Nat. Phys. 8, 422–428 (2012).
 42.
Cao, Y., Wang, H., Ouyang, Q. & Tu, Y. The freeenergy cost of accurate biological oscillations. Nat. Phys. 11, 772–778 (2015).
 43.
Sartori, P. & Pigolotti, S. Kinetic versus energetic discrimination in biological copying. Phys. Rev. Lett. 110, 188101 (2013).
 44.
Wang, S.W., Lan, Y. & Tang, L.H. Energy dissipation in an adaptive molecular circuit. J. Stat. Mech. 2015, P07025 (2015).
 45.
Sekimoto, K. Stochastic Energetics, vol. 799 (Berlin Springer Verlag, 2010).
 46.
Seifert, U. Stochastic thermodynamics, fluctuation theorems and molecular machines. Rep. Prog. Phys. 75, 126001 (2012).
 47.
Kubo, R. The fluctuationdissipation theorem. Rep. Prog. Phys. 29, 255 (1966).
 48.
Sethna, J. Statistical Mechanics: Entropy, Order Parameters, and Complexity, vol. 14 (Oxford University Press, 2006).
 49.
Matsuoka, K. Sustained oscillations generated by mutually inhibiting neurons with adaptation. Biol. Cybern. 52, 367–376 (1985).
 50.
Mello, B. A., Shaw, L. & Tu, Y. Effects of receptor interaction in bacterial chemotaxis. Biophys. J. 87, 1578–1595 (2004).
 51.
Tang, M. et al. Evolutionarily conserved coupling of adaptive and excitable networks mediates eukaryotic chemotaxis. Nat. Commun. 5, 5175 (2014).
 52.
Lindner, B., GarcíaOjalvo, J., Neiman, A. & SchimanskyGeier, L. Effects of noise in excitable systems. Phys. Rep. 392, 321–424 (2004).
 53.
Izhikevich, E. M. Dynamical Systems in Neuroscience (MIT Press, 2007)
 54.
Weber, A., Prokazov, Y., Zuschratter, W. & Hauser, M. J. Desynchronisation of glycolytic oscillations in yeast cell populations. PLoS ONE 7, e43276 (2012).
 55.
du Preez, F. B., van Niekerk, D. D., Kooi, B., Rohwer, J. M. & Snoep, J. L. From steadystate to synchronised yeast glycolytic oscillations i: model construction. FEBS J. 279, 2810–2822 (2012).
 56.
Gustavsson, A.K., Adiels, C. B., Mehlig, B. & Goksör, M. Entrainment of heterogeneous glycolytic oscillations in single cells. Sci. Rep. 5, 9404 (2015).
 57.
Chen, C., Liu, S., Shi, X.q, Chaté, H. & Wu, Y. Weak synchronisation and largescale collective oscillation in dense bacterial suspensions. Nature 542, 210–214 (2017).
 58.
Kawaguchi, K., Kageyama, R. & Sano, M. Topological defects control collective dynamics in neural progenitor cell cultures. Nature 545, 327 (2017).
 59.
Devreotes, P. N. et al. Excitable signal transduction networks in directed cell migration. Annu. Rev. Cell Dev. Biol. 33, 103–125 (2017).
 60.
Janetopoulos, C., Jin, T. & Devreotes, P. Receptormediated activation of heterotrimeric gproteins in living cells. Science 291, 2408–2411 (2001).
 61.
Loomis, W. F. Cell signaling during development of dictyostelium. Dev. Biol. 391, 1–16 (2014).
 62.
Danino, T., MondragónPalomino, O., Tsimring, L. & Hasty, J. A synchronised quorum of genetic clocks. Nature 463, 326–330 (2010).
 63.
Liu, C. et al. Sequential establishment of stripe patterns in an expanding cell population. Science 334, 238–241 (2011).
 64.
PotvinTrottier, L., Lord, N. D., Vinnicombe, G. & Paulsson, J. Synchronous longterm oscillations in a synthetic gene circuit. Nature 538, 514 (2016).
 65.
Toda, S., Blauch, L. R., Tang, S. K., Morsut, L. & Lim, W. A. Programming selforganizing multicellular structures with synthetic cellcell signaling. Science 361, 156–162 (2018).
Acknowledgements
The authors thank Allon Klein and Kyogo Kawaguchi for helpful discussions and suggestions on the manuscript. The work is supported in part by the NSFC under Grant Nos. U1430237, 11635002 and U1530401, and by the Research Grants Council of the Hong Kong Special Administrative Region (HKSAR) under Grant No. 12301514 and C201415G.
Author information
Affiliations
Contributions
S.W.W. and L.H.T. designed research, performed research, analysed data, and wrote the paper.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Wang, SW., Tang, LH. Emergence of collective oscillations in adaptive cells. Nat Commun 10, 5613 (2019). https://doi.org/10.1038/s41467019135739
Received:
Accepted:
Published:
Further reading

Partial synchronisation of glycolytic oscillations in yeast cell populations
Scientific Reports (2020)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.