To track the signalling state of individual cells during morphogenesis, we imaged cells expressing the intensiometric cAMP signalling reporter Flamindo2 (Odaka et al., 2014; Hashimura et al., 2019) at an endogenous genetic locus. For simultaneous analysis of cell motion, the cells also expressed a nuclear marker (H2B-mCherry) for cell tracking (Corrigan and Chubb, 2014). The fluorescence of Flamindo2 has an inverse relationship to the intracellular [cAMP]. This allows visualisation of the onset and progression of cAMP signalling wave patterns over cell populations during early development (Figure 1A and Figure 1—video 1). To monitor the motion and signalling of individual cells in their population context, we tracked the central position of the nuclei and measured the internal [cAMP] from the mean Flamindo2 intensity of the cell body (Figure 1B, Figure 1—figure supplement 1A and Figure 1—video 2). Analysing all cell tracks simultaneously provided the means to study the onset and adaptation of signalling behaviour and collective motion of single cells in their population context (Figure 1C–E).

Figure 1 with 6 supplements see all

Download asset Open asset

Early during development, cells predominantly fired asynchronously. The proportion of spontaneously activating cells declined as the signalling behaviour became more correlated, with periodic waves of signalling including an increasing proportion of the cells (Figure 1D). Prior to the onset of cAMP waves (around 4 hr after starvation), cells produced cAMP pulses spontaneously at a mean rate of around 1 pulse per 12 min, with a minimum of around 1 pulse per 4 min (Figure 1E). During this time, the population distribution of inter-activation intervals follows a gamma distribution (Figure 1E and Figure 1—figure supplement 1B). We interpret this fit as there being a minimum waiting time between pulses (a refractory phase), with the exponential tail of the distribution suggesting the probability of firing is independent of the previous firing event.

The initial cAMP waves (around 4–5 hr) are broad bands of broken waves amalgamating from numerous circular patterns that expand from spontaneously activating cells (Figure 1A). Activation events prior to synchronisation are spontaneous, spatially variable, and mostly do not initiate a cAMP wave. Together with the independence of successive firing events, these observations are inconsistent with the model that individual ‘pacemaker’ cells act as the organising source of global cAMP waves. However, it may be possible that cell clusters (arising for example from random motion) may collectively act as a pacemaker (Mirollo and Strogatz, 1990; Gregor et al., 2010). For example, if the cell cluster is sufficiently dense such that at least one cell in the cluster is likely to activate just after the refractory time. Nevertheless, any given cell will experience circular waves from several different sources with a variable periodicity (5–9 min).

Circular waves produce cAMP signals at irregular intervals which are typically greater than the inter-pules interval (Figure 1D). Consequently, not all cells will relay a circular wave because cells that fire synchronously with one wave can potentially reactivate and so be refractory for the next wave. This is highlighted by a multi-peaked population distribution for inter-activation intervals where the smallest peak is associated with spontaneous cAMP pulses, while the other peaks are associated with synchronised firing during a wave (Figure 1E and Figure 1—figure supplement 1B). Rotating spiral cAMP waves form from broken waves around 5–6 hr following starvation (Figure 1A). To interrogate how waves break and give rise to spiral waves, we engineered cells to express blue-photoactivable adenylyl cyclase (bPAC; Stierl et al., 2011) to control the time and place of cAMP waves by exposing a subset of cells to blue or green light (Figure 1—figure supplement 2 and Figure 1—video 3). In the absence of light pulses, we observed a steady train of periodic waves (Figure 1—figure supplement 2B). Photoactivation of a small patch of cells generated a circular wave which fused with oncoming waves (Figure 1—figure supplement 2A first and second row and Figure 1—figure supplement 2C). Wave breaks (Figure 1—figure supplement 2A forth – sixth row) occurred only when the oncoming wave collided with an asymmetric wave generated by the activated patch (Figure 1—figure supplement 2A third row). We interpret the asymmetry arising when cells activate at the border of two regions; one where cells are in a refractory state (preventing wave propagation) and the other in which cells are recovered (allowing wave propagation; Figure 1—figure supplement 2D). This experiment suggests that spiral waves might naturally arise from spontaneous cell activation behind a circular wave, where neighbouring cells are transitioning from a refractory to rest state.

Rotating spiral patterns generate waves that are more frequent and regular than circular patterns (Figure 1 5-6 hr onwards) and hence become the dominant pattern over time, overriding the less frequent circular waves, as reported in previous studies (Lee et al., 1996). The mean wave frequency steadily increased towards the maximum firing frequency of cells, around once every 3–4 min (Figure 1E). Consequently, the proportion of spontaneously activating cells tends to zero as there is no time to activate between waves (Figure 1D), and the population distribution of inter-pulse intervals tends to a normal distribution (Figure 1E). To understand the basis of the heterogeneity in inter-pulse intervals, we quantified the single-cell differences in activation times relative to the wave centre (Figure 1—figure supplement 3A, B). Aided by Markov modelling, we observed a short-term bias in firing tendencies: for example, a cell that fires ahead of its neighbours has a 60% chance of repeating this early firing, 35% chance of switching to late firing, and a 5% chance of skipping the subsequent wave (Figure 1—figure supplement 3C, D). Whether a cell fires before or after the wave does not affect its probability of skipping a wave. Overall, these data suggest no intrinsic difference between cells beyond a slight history dependence.

To relate how information encoded by spiral waves coordinates collective motion, we compared changes in cAMP wave dynamics to cell chemotaxis (Figure 2 and Figure 1—video 2). Chemotaxis is directed by the external cAMP concentration field, which is determined by the release of internal cAMP (reported by Flamindo2) from individual cells. The spatial organisation and motility of cells changed rapidly, occurring around 7 hr following starvation, coincident with the dominance of spiral waves (Figure 2A and Figure 2—figure supplement 1). During this process, cAMP waves steadily became more frequent, slower and thinner, and cells became faster, more directional, and persistent in their movement (Figure 2B and C). The profile of the mean cAMP pulse in single cells remained constant with each wave. Since the duration of a cAMP pulse ($\tau$) remains constant, and the wave speed ($v$) decreases, then the wave width ($w=v\times \tau$) also decreases. A thinner wave generates a steeper external cAMP gradient, which may enhance chemotaxis (Song et al., 2006).

Figure 2 with 1 supplement see all

Download asset Open asset

To quantify collective cell motion, we used the following two measures: (i) the mean chemotaxis index, the proportion of cell movement in the direction of the oncoming wave (Skoge et al., 2014; Figure 2B and D), and (ii) the local correlation of cell velocity (De Palo et al., 2017; Figure 2E). Cells consistently moved towards the wave source at wave fronts but remained spatially uniform before 6 hr of starvation due to random cell motion at the wave rear (Figure 2B bottom left panel). This is highlighted by a peak in the local correlation of cell motion (Figure 2E) and the mean chemotaxis index (Figure 2D bottom panel) at the wave front, both of which then decay towards zero (indicating random cell motion). For infrequent, fast and wide waves (i.e. circular patterns), the chemotaxis index and local correlation of cell motion reached their baseline values at the following wave front, indicating cells move randomly. This random cell motion disrupts the order in collective cell motion initially induced by the front of the wave. For frequent, slow and thin waves (i.e. rotating spiral waves), the chemotaxis index and velocity correlation significantly change as the cAMP wave frequency increases between 6 and 8 hr following starvation. From 6.5 hr onwards, when the wave period is less than 4.5 min, the measures of collective cell motion do not relax back to their previous values by the next wave and hence gradually increase over time to constant, elevated values across the entire wave (Figure 2D and E). In parallel, the cell speed at the front of the wave increases incrementally with each wave until it peaks at around 7 µm min⁻¹ (Figure 2D second panel). These data are consistent with models where decreasing the wave period causes cells to gradually become more persistent towards the direction of the wave source and hence gain consistent collective motion (Skoge et al., 2014).

To test the relationship between wave frequency and collective migration in the full population context, we photoactivated bPAC in a small patch of cells at different intervals and monitored the resultant activation waves (Figure 3 and Figure 3—video 1). Consistent with our observations shown in Figure 2, cells collectively migrated towards the signalling centre only when the induced wave period was less than 4.5 min (Figure 3A). However, photoactivation of cells with a period less than 3.5 min did not result in aggregation. We infer that this is because cells cannot relay a steady periodic signal close to their maximum firing rate initially, resulting in a missed wave every second photoactivation event (Figure 3C), leading to the signalling centre becoming ‘taken over’ by another set of period waves. Motivated by the observed change in wave frequency during aggregation, we sought to test if it were possible to induce large cell aggregates by photoactivating cells at a faster rate than what we observe experimentally by reducing the wave period gradually from 5 min to 3.5 min from 5 to 6 hr following starvation, rather than 6–8 hr (Figure 3C second panel). This was indeed achievable, resulting in the accumulation of a significantly larger number of cells around the signalling centre and conglomerates of several aggregates (Figure 3B). These artificial signalling centres induce large aggregates by entraining the territories of neighbouring signalling centres. Our data support the notion that the increase in wave frequency is a crucial determinant of the transition to collective motion. We now turn to the question of how cells generate this increasing wave frequency.

Figure 3 with 1 supplement see all

Download asset Open asset

To characterise the regulation of spiral wave dynamics, we imaged and analysed the motion and signalling of individual cells around the spiral core following the onset of spiral formation. The spiral core is the region about which the spiral tip (the end of the activation wave) moves. Based on the general theory of excitable media, the spiral tip will sustain itself and the entire spiral wave by virtue of circular motion (Tyson and Keener, 1988; Winfree, 2001). In other words, the circulation of the spiral tip is the organising source of the spiral wave. To better understand the dynamics of spiral tip circulation, and how it relates to the changes in wave dynamics, we tracked the position of the spiral tip and the shape of the spiral wave from spiral formation to aggregation (Figure 4, Figure 4—figure supplement 1, Figure 4—videos 1–3). To track tip motion, we mapped the [cAMP] of individual cells to a phase variable measuring the sequence of changes in the internal [cAMP] during a cAMP pulse (Figure 4B, Figure 4—figure supplement 1B, C). During a cAMP pulse, this variable increases from 0 (baseline [cAMP]) to ±π (maximum [cAMP]) and then back to 0. This approach reveals the spiral singularity — or topological defect — which is the point in space around which cells are distributed across the full spectrum of phase states. At each time point, we estimated the spiral tip location as the mean position of maximally activated cells closest to the singularity (Figure 4B). Additionally, we fitted a spiral curve to the spatial distribution of maximally activated cells (Figure 4A and C first panel).

Figure 4 with 4 supplements see all

Download asset Open asset

Spiral waves significantly change in rotation rate and curvature during aggregation, which we refer to as spiral wave progression (Figure 4C and Figure 4—videos 1–3). These changes in spiral wave structure and dynamics are coupled to changes in spiral tip circulation (Figure 4B–D). Each rotation of the spiral tip generates a wave: the rate at which the activation waves circulate around the signalling centre sets the frequency of global activation waves disseminated from the signalling centre. The spiral tip rotational period decreased from 5.5 min to 3.5 min. These changes in the rotational period of the spiral tip translate to the changes in the rotation rate of the spiral wave: the decrease in the spiral tip rotation period (Figure 4D first panel) is consistent with the decrease in the wave period away from the spiral core that coincides with the onset of collective migration (Figures 2C and 3). The acceleration of spiral tip circulation is associated with spiral core contraction; the radius of circulation steadily reduced from around 100 µm to 20 µm (Figure 4D second panel). Meanwhile, during this contraction, the spiral tip speed decreased from 100 to 120 to 40–50 µm min⁻¹ (Figure 4D third panel).

How do cells at the signalling centre drive spiral wave progression? To measure the fundamental signalling parameters and how they change over time, we utilised the general theory of excitable media (Tyson and Keener, 1988). Dictyostelium cAMP waves can be shown to satisfy this theory (Tyson et al., 1989; Foerster et al., 1990). Our high-resolution datasets now make it possible to apply this theory to understand how spiral waves change during aggregation, independently of the specific equations that model cAMP release, which are still disputed (Sgro et al., 2015; Kamino et al., 2017). According to this theory, the shape of the spiral wave must satisfy a curve that is approximately Archimedean (Materials and methods, Equation 2). The spiral curve depends on a single parameter $k=N_0/\epsilon D$ (referred to here as the curve parameter) which relates the signalling parameters $N_0$ (the planar wave speed), $D$ (the signal diffusion constant), and $\epsilon$ (the ratio of the rates of cell inactivation and activation, which we call the cell activation parameter; Materials and methods, Equation 3). Crucially, the inverse of $\epsilon$ is the measure of cell excitability and specifies the rate at which cells activate. By fitting the spiral curve to the spatial distribution of maximally activated cells at each time point (Figure 4A; Figure 5A, Figure 4—videos 1–3), throughout spiral wave progression, we obtain a measurement of how the cell excitability and planar wave speed change over time (Figure 5B second column). Specifically, both measures of the cell signalling phenotype can be expressed in terms of spiral wave properties that we can measure (Materials and methods, Equation 4), such as the normal wave speed and wave curvature (Figure 5B first column).

Figure 5 with 1 supplement see all

Download asset Open asset

This analysis reveals that the activation parameter and planar wave speed both decrease during spiral wave progression (Figure 5B). Firstly, we find that the planar wave speed halves in value as the wave period decreases towards its minimum (Figure 5D). The theory of excitable media provides the following interpretation: as the wave period decreases, cells have less time to reset between waves, causing slower cell activation and wave speed (the dispersion relation), which gives rise to the deceleration of wave speed during core contraction. This deceleration quantitatively matches the changing wave dynamics away from the spiral core (Figures 2C and 5B), meaning the circulation of the spiral tip determines the wave speed away from the core, in addition to the frequency. Secondly, the theory interprets the increase in spiral wave rotation rate and curvature to be a consequence of cells becoming more excitable over time (Figure 5C). In other words, the spiral waves progress because cells can activate faster. Indeed, although there was temporal variation between experimental replicates, the activation parameter collapses onto a highly consistent curve when plotted against the wave period (Figure 5C), showing how the wave dynamics are strongly coupled to the excitability of the cells. A simple quadratic provides a good fit regardless of derivation used for quantifying the activation parameter (Materials and methods, Equations 4 and 5). In summary, this analysis implies the change in cell excitability is a key element driving spiral wave progression. We now evaluate different hypotheses for mechanisms driving spiral wave progression.

We have investigated two hypotheses for spiral wave progression. The prevailing paradigm (Hypothesis 1) is that excitability increases due to an increasing cell density; however, this model is not supported by our data comparing density with wave period and excitability (Figure 6). Our reformulation of this hypothesis predicted that cell rearrangements drive the increase in the spiral wave rotation rate via contraction of spiral tip circulation (Figure 7). In contrast, Hypothesis 2 proposes that the wave dynamics change cell excitability (Figure 8). The key difference between these hypotheses is that in Hypothesis 1, spiral tip contraction is the driver of spiral progression (Figure 7F), whereas in Hypothesis 2, contraction is a consequence (Figure 8—figure supplement 1). Exploring the properties of the models underlying these hypotheses provides a clear experimental test: to constrain the radius of spiral tip circulation.

To experimentally test the competing hypotheses, we imaged Dictyostelium aggregation in the presence of physical blocks (low melt agarose particles with radii ranging between 25 and 150 µm; Figure 9 and Figure 9—video 1). Spiral waves can become pinned to the blocks (Figure 9A), with their minimum radius of circulation set by the block radius (Figure 9B). Spiral waves that are pinned to a block cannot fully contract (Figure 9C first panel) but nevertheless undergo a steady increase in rotation rate and hence decrease in wave rotational period (Figure 9C third panel). The circulation rate increases because the spiral tip speed does not decrease to the same extent as the unconstrained spirals (Figure 9C second panel). This observation is in agreement with the model simulations shown in Figure 8, revealing the same effects as in the data – wave speed increased to compensate for larger radius of circulation, generating a reduced rotational period. This observation also reflects our simulations predicting that spiral waves can progress in the presence of physical blocks due to a coupling between the cell excitability and wave dynamics. In addition, this coupling between wave frequency and excitability is also apparent in our photoactivation experiments, in which a gradually increasing pulse frequency enhanced aggregate formation (Figure 3). Together, these observations are consistent with Hypothesis 2, but not Hypothesis 1: the contraction of the spiral tip circulation is a consequence, not a driver, of spiral wave progression.

Figure 9 with 1 supplement see all

Download asset Open asset

Overall, our approach provides a clear interpretation of our data that spiral wave progression results from positive feedback between cell excitability and wave frequency, arising only in the context of the spiral wave pattern (Figure 8—figure supplement 1). That is, from the initial symmetry breaking event (Figure 1—figure supplement 2), the formation of a rotating spiral wave produces stable periodic activation waves that increase cell excitability, then increases the spiral rotation rate and hence wave frequency, which then feedback onto further increases in cell excitability.

To conclude, we now describe the cAMP wave dynamics following aggregation. As development proceeds, more and more cells move into the aggregate. Previous studies show the activation waves are organised as a multi-armed spiral in aggregates (Siegert and Weijer, 1995). How does this continued cell accumulation affect the cAMP dynamics at the signalling core and, consequently, the encoding of information by cAMP waves?

To characterise the single- to multi-armed spiral transition, we monitored the structure and dynamics of the activation waves following ring collapse (Figure 10, Figure 10—figure supplement 1 and Figure 10—video 1). Single-armed spirals destabilise only once the cellular ring fully contracts. Contraction results in a phase singularity (convergence of all signalling phases) that is fixed at the aggregate centre (Figure 10—figure supplement 1A). The phase singularity and spiral wave are unstable and disintegrate over time. This is likely because the cells are not infinitely small and can only be in one signalling phase at a time. In other words, the cells at the centre cannot activate at a rate necessary to support the full spectrum of signalling phases. The destabilisation process begins by cells signalling out-of-phase with the spiral wave, resulting in the emergence of transient patches of cell activation that gradually increase in number, size, and lifetime (Figure 10, around 2 hr following spiral formation). These activation patches rotate about the aggregate while fusing with, and breaking apart from, other activation patches and the original spiral wave, which progressively disintegrates (Figures 10A and 2–3 hr following spiral formation). This process concludes with the appearance of three stable, slow-moving and uncurved activation waves that rotate about the same fixed point at the aggregate centre (Figure 10A, from 3 hr following spiral formation). At this point, the signalling regime can be classed as a multi-armed spiral. Although the ring has collapsed, the centre of the aggregate contains relatively few cells (Figure 10—figure supplement 2). At this point, there is no phase singularity: the few cells that occupy this central region appear to have signalling phases that are uncorrelated with each other and the surrounding multi-armed spiral wave (Figure 10—figure supplement 1 – last panel). This central region of uncorrelated signalling phases is characteristic of theoretical chimera spiral waves (Totz et al., 2018).

Figure 10 with 3 supplements see all

Download asset Open asset

To analyse the transition from the single- to multi-armed spiral, we generated a set of kymographs of cell activation states across the aggregate: one that is cross-sectional and the others that are circular (Figure 10B). The cross-sectional kymograph shows a disordered signalling regime between 2 and 3 hr following spiral formation, from around 45 min after cellular ring closure. In the circular kymographs, both the single- and multi-armed spiral waves appear as diagonal lines across space and time (<2 and >3 hr), meaning a constant rotation rate. These data show that the single-armed spiral (<2 hr) has a rotational period of 3.4±0.2 min (the minimum cell firing period, see Figures 1, 3 and 4), which is three times faster than the three-armed spiral whose rotational period is 10.1±0.7 min and wave period is 3.3±0.5 min. This means that the overall frequency with which a cell experiences a wave is the same because the three-armed spiral provides three waves in its longer rotation period. As such, the frequency of cAMP waves disseminated from the signalling centre, along the streams, remains the same for both signalling patterns. That is, both patterns encode the same temporal information.

Slow- and fast-waves coexist at the outside of the aggregate midway through the transition to the multi-arm spiral. This is visualised as a mixture of lines of two different gradients in the kymograph (Figure 10B, fourth panel, 2–3 hr). This observation suggests that both signalling regimes coexist by the mixing of cells that relay either the single- or multi-armed spiral wave. The multi-armed spiral wave eventually becomes dominant as the proportion of cells that relay the single-armed spiral wave declines. Altogether, our data suggest that the signalling pattern adapts to the change in aggregate geometry while maintaining an active phase wave that generates waves at the maximum frequency.

Dictyostelium AX3 cells were engineered to express Flamindo2 and H2B-mCherry from safe-haven genomic locations for stable and uniform expression. Flamindo2 is a cAMP biosensor (Odaka et al., 2014; Hashimura et al., 2019) with a fluorescence intensity that is inversely proportional to internal cAMP levels. A codon-optimised sequence encoding the Flamindo2 reporter was knocked-in to the act5 locus for stable and uniform expression (Paschke et al., 2018; Tunnacliffe et al., 2018). H2B-mCherry, knocked into the rps30 locus, labels the nucleus for cell tracking at highcell densities (Corrigan and Chubb, 2014). For optogenetic activation of cAMP signalling, we used the light-activatable adenylyl cyclase, bPAC (Stierl et al., 2011). A codon-optimised bPAC cDNA was cloned into the extrachromosomal expression vector, pDM1203, then transfected into the act5-Flamindo2 cells.

Cells were grown at 22°C in 12 ml of Formedium HL5 medium supplemented with penicillin and streptomycin in 10 cm tissue culture dishes. Cultures were diluted to 25% of their density each day to stop cell entry into stationary phase. Cells were not used beyond 9 days of culture. For development, log-phase cells at 80% confluency were washed and resuspended in KK2 (20 mM KPO₄ pH 6.0) and then plated on 1.5% agar at between 4×10³ and 5×10³ cells mm⁻². After 1 hr in a humidified chamber, a 1.2 cm × 1.2 cm slab of agar was cut and gently placed cell-face-down onto a glass imaging dish (Bioptechs deltaT) and then immediately submerged in silicon oil to prevent desiccation. To incorporate physical blocks into the cell aggregation field, we sprinkled low-melting point agarose powder onto the KK2 agar bed before the agar had set.

Cells were imaged at 22°C on a custom-built inverted wide field microscope (Cairn Research) with 470 nm and 572 nm LED light sources (Miermont et al., 2019). This was equipped with a Prime 95B CMOS camera (Photometrics) and 10× UplanFL N objective (Olympus). A 1.2 mm × 1.2 mm area was imaged at resolution 1 µm² per pixel every 4 s for 8 hr from 2 hr after starvation. Frame stitching was used to expand the effective field-of-view size. To image cell motion and signalling about the spiral core, including the experiments with the physical block, we imaged a 1.2 mm × 1.2 mm area for 2 hr following from 5 to 6 hr after starvation. For this, we monitored cAMP wave patterns in real time and identified waves that had started to rotate about an axis.

For light activation of cAMP signalling, we used bPAC-expressing act5-Flamindo2 cells and a 3i spinning-disk confocal microscope with a 10× objective and Prime 95B CMOS camera (Photometrics). To initiate periodic cAMP waves, cells contained within a circular area with a ≈100 µm diameter were pulsed with 488 nm light for 1 s, at intervals to mimic physiological cAMP pulsing. Flamindo2 images (excitation 515 nm) were collected throughout the pulse sequence, every 5 s. To initiate a rotating spiral wave, we repeated the photoactivation protocol above, except that cells were photoactivated at operator-controlled points in time and space. We monitored waves in real time and identified periodic waves away from the signalling centre (around one per 5 min, 5 hr following starvation). To induce circular wave patterns, we photoactivated cells around 3 min after the previous wave. To induce asymmetric wave patterns (resulting in the broken and hence spiral wave), we photoactivated around 2 min after the previous wave.

The following image analysis pipeline was implemented in Matlab 2020b. To track the position of individual cells from image sequences, the IDL particle tracking algorithm (Crocker and Grier, 1996) was applied to the cell nuclei channel (H2B-mCherry) of the image time series. The velocity of each cell at each time point was estimated by the first derivative, second order, central finite difference of cell position. To measure the mean fluorescence intensity of each tracked cell, watershed segmentation (the Fernand Meyer algorithm) seeded from the tracked cell positions was applied to the cell body channel (Flamindo2). To infer internal cAMP levels (non-dimensional), the tracked mean fluorescent intensity of each cell was normalised and smoothed using Matlab’s Savitzky-Golay finite impulse response smoothing filter. Two measures were used to identify cAMP pulses: (i) the similarity between the normalised intensity and a square pulse of length 2 min and (ii) the signalling phase $\phi$. Variable $\phi$ was estimated by stratifying the time series of the internal cAMP level $c$ by its time derivate $\dot{c}$ (estimated by the first derivative, second order, central finite difference). The value of $\phi$ is set as the angle around the point $c=c^{*}$ and $\dot{c}=0$:

In addition to the time series of the inferred cAMP levels, we recorded the time points of the maximal cAMP level during each pulse, and the time points of cell activation (${\displaystyle \dot{c}>0}$) and deactivation (${\displaystyle \dot{c}<0}$), measured by the time point where the cAMP levels are half-maximal.

To quantify the internal [cAMP] signalling and phase field following aggregation, we tracked the fluorescence intensity of each individual pixel as a function of time, over the entire 2D image sequence. The fluorescence intensity time series was oscillatory and typically increased over time due to cell crowding. We used Matlab’s Savitzky-Golay finite impulse response smoothing filter to separate the oscillatory and background component of the signal. The background component was used to estimate the cell density of the aggregate. The oscillatory component was used to estimate the spatial distribution of signalling states and phases across the aggregate, using the same technique described above.

To quantify the evolving dynamics of the spiral wave from cell track data, we first estimated the location of the spiral phase singularity at each time point as the central point of a circular region of radius 100 µm with the largest variance in internal [cAMP] across each cell. The location of the spiral tip was estimated by the mean position of five maximally activated cells closest to the location of the phase singularity. We estimate the axis of rotation at each time point by fitting a circle to the track of the spiral tip, half a revolution backwards and forwards in time. Using the Matlab function ‘fit’, we fit the spiral curve given by Equation 2 (see below) to the spatial distribution of maximally activated cells, starting from the axis of spiral tip circulation (Tyson and Keener, 1988). To achieve this, we determined the magnitude ($r$ in Equation 2) and orientation ($\theta$ in Equation 2) of the vector extending from the axis of spiral rotation to the location of each maximally active cell. This fit provides the spiral orientation (${\theta }_0$ in Equation 2) and curve parameter ($k$ in Equation 2).

Away from the spiral core, the time and place of activation for each cell were clustered in three dimensions (two spatial and one temporal) using the DBSCAN (density-based clustering of applications with noise) algorithm as a means to identify cells with correlated signalling states (e.g. during the propagation of cAMP waves) and also, by exclusion, the cells that activate independently. This approach reveals the clusters of cell activation associated with each cAMP wave. For each wave, we estimated the dynamic curve of the wave by fitting an equation for a cone ($t-t_0=Nr,r=\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}\mathrm{}$, where $t_0$ and ${(x}_0,y_0)$ is spatiotemporal origin of the wave) to the time evolution of the spatial distribution of maximally activated cells. The slope of each of the cones fitted to the spatial distribution of maximally activated cells at each time point was used to deduce the speed of each wave.

To estimate the mean cell velocity field and density around the spiral wave, we measured the mean velocity and total number of cells at regular lattice areas (step size 10 µm) in the reference frame that rotates with the fitted spiral parameter ${\theta }_0$ (Figure 6—video 2). The streamlines across this field were analysed using the Matlab function ‘streamlines’.

To analyse the behaviour of single cells during stable periodic waves, the fitted cone equations were used to estimate the mean internal [cAMP] and velocity of each cell as a function of distance from the wave centre. This allows an estimation of measures such as the angle difference between the direction of wave travel and cell motion, for each cell:

where $\widehat{{\mathit{v}}_i}$ is the normal velocity of cell $i=1,\dots ,N$, and $\widehat{\mathit{w}}\left({\mathit{X}}_i\right)$ is the normal velocity of the cAMP wave at the location of cell $i$, ${\mathit{X}}_i$ . The following two formulae were used to measure collective cell motion: (i) the mean chemotaxis index (the proportion of cell movement in the direction of the oncoming wave) as a function of distance from the wave:

where ${\Omega }_j$ represents the subset of cells (subset size $N_{{\Omega }_j}$) located at a distance $j\times 10$ µm ($j=\dots -\mathrm{1,0},1,\dots$) from the wave centre, and (ii) the local correlation of cell velocity:

where ${\Phi }_i$ denotes the subset cells (subset size $N_{{\Phi }_i}$) within a 50 µm radius of cell $i$ (De Palo et al., 2017). Local cell density was measured by the mean density of cells within a 60 µm radius of each cell.

This section outlines the elements of the theory we have used, and we used these to interpret spiral wave progression (Tyson and Keener, 1988; Winfree, 2001).

The eikonal relation dictates how the normal speed along the wave front $N$ changes due to wave curvature $K$:

where $N_0$ is the planar (no curvature) wave speed, $D$ is the signal diffusion coefficient (which we assume remains constant), and $\epsilon$ is the ratio of the rates of cell inactivation and activation (which we call the cell activation parameter; Keener, 1986). For a wave extending to the pivotal axis, the solution to Equation 1 is accurately approximated by the following curve (in polar coordinates, $x=r\cos \theta (r,t),y=r\sin \theta (r,t)$):

where $m$ is a constant. Parameter $k$ (the curve parameter) defines the curvature of the wave, and ${\theta }_0$ is the spiral angle of rotation, with rotation rate, ${\omega =\dot{\theta }}_0$ (Tyson and Keener, 1988). We fit Equation 2 to the spatial distribution of maximally activated cells at each time point, throughout the entire progression of the spiral wave. The fit is accurate for all time points, which validates the use of the theory of excitable media even though cells are discretely distributed across space. This fit provides a measurement of the time evolution of the curve parameter $k$, rotation rate $\omega$, and hence the wave period $T=2\pi /\omega$ (Figure 4D). The value of $k$ relates the planar wave speed to the cell excitability:

Substituting Equation 3 into Equation 1 for $N_0$ allows the cell activation parameter $\epsilon$ and planar wave speed $N_0$ to be expressed in terms of variables that can be directly measured from data:

Consistent with the traditional paradigm of spiral wave progression (Höfer et al., 1995; van Oss et al., 1996), we find that the cell activation parameter, which is the inverse of cell excitability, decreases during spiral wave progression (Figure 5). This result illustrates how we can identify key signalling properties of the cells from their population signalling patterns, at each time point during aggregation. Plotting $N_0$ from Equation 4 against the wave period $T$ provides a prediction of the dispersion relation, which quantifies how a decrease in wave period reduces the wave speed because cells do not fully reset between waves, causing them to activate slower. The curvature relation describes the relationship between radius and rate of circulation of the spiral tip:

where $r_0$ is the radius of spiral tip circulation, which can measure from our data (Tyson and Keener, 1988). In support of our approach, the curvature and eikonal relations produce similar values of the cell activation parameter at all time points. We use the curvature relation to predict the circulation radius of the spiral tip based on our spiral measurements.

A hybrid PDE ABM was developed to study the interplay between cell motion, spiral wave dynamics, and spatial patterning of cells around the spiral core.

External cAMP is secreted by active cells. We define active cells as having a positive activator value ${\displaystyle A_i>0}$ (Sgro et al., 2015). Once secreted, we assume external cAMP diffuses and decays. The 2D external [cAMP] field $c$ was modelled by the following reaction-diffusion equation:

where $\beta$ and $\sigma$ are the rates of external cAMP decay and release from individual cells, $D$ is the external cAMP diffusion coefficient, $N$ is the number of cells, $H$ is the Heaviside function, and $\delta$ is the Dirac delta function, modelling the point source release of cAMP from cell $i$ (${\displaystyle 1\le i\le N}$), with position ${\mathit{X}}_i$ and activation state $A_i$. The spatial domain is $\mathrm{\Omega }={\left[-\frac{L}{2},\frac{L}{2}\right]}^2\subset {\mathbb{R}}^2$ where ${\displaystyle L>0}$ is the domain size. We used same parameter values of Vidal-Henriquez and Gholami (Vidal-Henriquez and Gholami, 2019). The following modified FitzHugh-Nagumo model was used to describe the excitable release of cAMP from each cell:

where $I_i$ is the inhibitor level of cell $i$, and the value of each parameter was the same as Sgro et al., 2015, with exception of the time scale $\tau$ which was adjusted to match experimental data.

We tested several mathematical models of cell chemotaxis against the data for the velocity of individual cells as a function of the external cAMP field. The external cAMP field was inferred in our experimental data by substituting the measured position and discretised activation state of each cell in the reaction-diffusion equation above, which provides a prediction for the external [cAMP] at the location of each cell $c_i=c_i\left({\mathit{X}}_i\right(t),t)$ and the spatial gradient $\left|c_i\right|$ (Figure 7—videos 1 and 2). The model that best fit the experimental data was a simple system of linear ordinary differential equations that encapsulate chemotactic memory:

where ${\displaystyle 0\le {\theta }_i<2\pi }$ and ${\displaystyle v_i\ge 0}$ are the direction and speed of cell movement, and ${\displaystyle s_i\ge 0}$ is the sensation variable to external cAMP as used in Höfer et al., 1995. This variable causes cell chemotaxis to be sensitive to the external cAMP gradient at the front of the wave, but not the rear. Variable $\left|\nabla c_i\right|$ and ${\displaystyle 0\le {\theta }_i^c<2\pi }$ are, respectively, the magnitude and direction of the external [cAMP] spatial gradient evaluated at ${\mathit{X}}_i$. Model parameters ${\alpha }_m$, ${\beta }_m$, ${\gamma }_m$, ${\sigma }_m$, and ${\lambda }_m$ were determined during the model identification using the Matlab function ‘fminsearch’ to obtain the best fit to experimental data. The first equation models the turning of the cell towards the local gradient of external cAMP. The turning speed depends on a cell’s sensitivity and the size of the cAMP gradient. The second equation describes adaptation of the cell speed. It drives cell speed to be proportional to $s_i\left|\nabla c_i\right|$ , i.e., higher sensitivity and larger local gradients will lead to faster speeds. Finally, the third equation models sensitivity: If a cell experiences a high concentration of external cAMP $c_i$, it will become de-sensitised, and it takes some time to become sensitive again. The variable $s_i$ provides ‘inertia’ or ‘memory’ that makes cell chemotaxis to be dependent on the history of external [cAMP] dynamics. When comparing the model predictions to the data, this effect was crucial to ensure that cells move towards the wave source at the wave front while also neglecting the wave rear. Specifically, between waves, the value of $s_i$ recovers to a state where the cells are sensitive to the external cAMP gradient at the front of the wave but subsequently $s_i$ reduces to a state where the cells are insensitive to the external cAMP gradient at the rear of the wave. For the final cell movement model, we also consider size exclusion effects between cells by including a cell-cell repulsion term. This means ${\theta }_i$ and $v_i$ are understood to mean the orientation and speed a cell would have if not in contact with other cells. We define as $d_{ij}=|{\mathit{X}}_i-{\mathit{X}}_j|$ the distance between the cell centres of the $i$-th and $j$-th cell. A cell’s movement now follows:

i.e., cells move with speed $v_i$ in direction ${\theta }_i$ in the absence of cell-cell repulsion. The second term describes a cell-cell repulsion with maximal repulsion magnitude ${\rho }_R$ where only cells within a distance of $d_R$ of each other exert a pushing force. See Table 1 for parameter values.

We use no-flux boundary conditions for the external cAMP concentration, i.e., cAMP cannot leave the domain. Similarly, cells cannot leave the domain. Since the dynamics cause a concentration of cells towards the spiral centre, cells are being depleted from the very edge of the simulation domain, and hence, we only considered the middle part of the simulation domain to represent the biological situation. For all numerical experiments, we initialise cell positions ${\mathit{X}}_i$ and orientations ${\theta }_i$ using a uniform random distribution on $\mathrm{\Omega }$ and $[0,\mathrm{}2\mathrm{\pi })$, respectively. Cell sensitivities $s_i$ were initially set to 1 and cell speeds $v_i$ to 0. The initial external cAMP concentration, activator, and inhibitor values were initialised as shown in Figure 7—video 3. The initial state was created by disabling cell movement and: (1) letting a vertical cAMP wave moving from left to right across a periodic domain until it is equilibrated, (2) setting the external cAMP in the upper-half of the domain to zero and setting the inhibitor value in the upper half on the domain to a high value (we used 2.5), and (3) setting the boundary conditions to no-flux and restarting the simulation.

Model equations were solved numerically using Matlab. The external [cAMP] field was solved across a 50 grid points per mm with the finite difference, implicit Euler method in time and central difference in space. Cell velocities and activation states were solved using the Matlab function ‘ode45’.

To test whether the coupling of cell excitability to the external [cAMP] dynamics sufficiently describes spiral wave progression, we extended a mathematical model of cAMP spiral waves to simulate the time evolution of the external [cAMP] field over time. We adopt the model by Tyson and Murray, 1989, which is of the general form of continuum excitable media with the dynamics of cAMP release governed by the Martiel-Goldbeter model:

where $u$ and $v$ (non-dimensional) are the spatial distribution of external [cAMP] (which can diffuse) and an internal inactivator (which cannot diffuse). We adopted the following non-dimensional parameter values from Tyson and Murray, 1989, most of which were rounded to the nearest order of magnitude: ${\alpha }_u={10}^3,{\beta }_u={10}^2,\kappa ={10}^1,L_1={10}^1,L_2={10}^{-1},{\lambda }_1={10}^{-4},{\lambda }_2=0.5,c={10}^1$. Equation 6 is dimensionalised by the time scale $\tau =0.075$ min and diffusion coefficient $D=3\times {10}^4$ µm² min⁻¹. We note that solutions to this equation rapidly tend to an equilibrium state. Motivated by our finding that the cell excitability is coupled to the cAMP wave dynamics, we extended Equation 6 by letting the cell activation variable, which we previously called the cell activation parameter (this inverse being a measure of cell excitability), to be a simple time-dependent function of the external [cAMP]:

Where ${\alpha }_{\epsilon }=\frac{1}{50},{\beta }_{\epsilon }=\frac{35}{50},{\epsilon }_0=10$. Equation 7 states that the cell activation variable tends to a baseline value ${\epsilon }_0$ and decreases in response to external cAMP. The simulations were performed on a 4.5 mm × 4.5 mm grid with 0.025 mm spacing and 0.05 min time step. The initial condition was the state where $=0,v=0$ , $\epsilon =\frac{{\epsilon }_0}{2}$, except for a rectangular region representing a broken wave, with $u=10$ and $v=0$ and to one side $u=0$ and $v=1$. To numerically solve Equations 6; 7, we used the implicit finite difference method for the diffusion term and the Runge-Kutta (Matlab ‘ode45’ function) method for the reaction terms. The Neumman boundary condition (no flux) was used.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication. For the purpose of Open Access, the authors have applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

This study was supported by a Wellcome Trust Senior Fellowship (202867/Z/16/Z) to JRC. The funder had no role in study design, data collection and interpretation, or the decision to submit the work for publication. We thank Elizabeth Westbrook, Ricardo Barrientos, Courtney Lancaster, Kees Weijer, and Allyson Sgro for comments on the manuscript.

© 2023, Ford et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.