Spontaneous center formation in Dictyostelium discoideum

TL;DR

This paper investigates how spontaneous formation of activity centers and target waves of cAMP in Dictyostelium discoideum can be explained by adding cell inhomogeneity to a reaction-diffusion model, revealing the role of cell density and wave dynamics.

Contribution

It introduces a mechanism for spontaneous center formation in D.d. by incorporating cell inhomogeneity into a reaction-diffusion model, advancing understanding of natural pattern emergence.

Findings

Stable activity centers emerge in high-density areas.

cAMP wave velocity increases with cell density.

Cells react to waves and form aggregation streams.

Abstract

Dictyostelium discoideum (D.d.) is a widely studied amoeba due to its capabilities of development, survival, and self-organization. During aggregation it produces and relays a chemical signal (cAMP) which shows spirals and target centers. Nevertheless, the natural emergence of these structures is still not well understood. We present a mechanism for creation of centers and target waves of cAMP in D.d. by adding cell inhomogeneity to a well known reaction-diffusion model of cAMP waves and we characterize its properties. We show how stable activity centers appear spontaneously in areas of higher cell density with the oscillation frequency of these centers depending on their density. The cAMP waves have the characteristic dispersion relation of trigger waves and a velocity which increases with cell density. Chemotactically competent cells react to these waves and create aggregation streams even with very simple movement rules. Finally we argue in favor of the existence of bounded phosphodiesterase to maintain the wave properties once small cell clusters appear.

Figures

Tables

Table 1. Table 1: Parameters used for simulations of equation ( 1 ) in the first row. Second row: Parameters used to approximate the production function in equation ( 2 ). Third row: Parameters used when cell movement was included

$c=10$	$h=5$	$k_1=0.09$ min-1
$\kappa =18.5$	$\sigma =0.55$ min-1	$k_i=1.7$ min-1
$k_t=0.9$ min-1	${ℒ}_1=10$	${ℒ}_2=0.005$
$q=4000$	${\lambda }_1={10}^{-4}$	${\lambda }_2=0.2575$
$\theta =0.01$	$D=0.024$ mm2/ min	$r=10$ $\mu$m
$a_0=0.04423$	$a_1=0.1526$	$a_2=328.2$
$v_c=0.02$ mm/min	${\rho }_c=0.6$	$\Delta {\gamma }_c=25.98$ nM/mm

Equations

$$k_1^{-1}{\partial }_t{\rho }_i$$

$${\partial }_t{\beta }_i$$

$${\partial }_t\gamma (x,y)$$

$$\centering\begin{array}[]{cccccc}\text{with}&f_{1}(\gamma)=\dfrac{1+\kappa\gamma}{1+\gamma},&f_{2}(\gamma)=\dfrac{\mathcal{L}_{1}+\kappa\mathcal{L}_{2}c\gamma}{1+c\gamma},&\Phi(\rho,\gamma)=\dfrac{\lambda_{1}+Y^{2}}{\lambda_{2}+Y^{2}},&\text{ and}&Y(\gamma,\rho)=\dfrac{\rho\gamma}{1+\gamma},\end{array}\@add@centering$$

$${\partial }_t\gamma =D{\partial }_{xx}\gamma -k_e\gamma ,$$

$$\gamma =A{e}^{-\sqrt{k_e/D}|x|}\text{if }|x|>L$$

$$0=a_0+(a_1-k_e)\gamma +a_2{\gamma }^2+D{\partial }_{xx}\gamma .$$

$${\partial }_{xx}\gamma =-\frac{\partial V}{\partial \gamma }\text{with}V=\frac{a_0}{D}\gamma +\frac{(a_1-k_e)}{2D}{\gamma }^2+\frac{a_2}{3D}{\gamma }^3$$

$${\int }_{\gamma }^{{\gamma }_M}\frac{d\gamma }{\sqrt{2[E({\gamma }_M)-V(\gamma )]}}=x.$$

$$\sqrt{\frac{k_e}{D}}\gamma (x=L)={\partial }_x\gamma {|}_{x=L}=\sqrt{2(E-V(x=L))},$$

$$k_1^{-1}{\partial }_t{\rho }_i$$

$${\partial }_t{\beta }_i$$

$${\partial }_t\gamma (x,y)$$

Full text

Introduction

Dictyostelium discoideum (D.d.) is a social amoeba that under adverse environmental conditions begins a developmental program in order to survive[1]. After 4-5 hours into this program, the cells start to produce the signaling chemical cAMP (cyclic adenosine monophosphate), which travels the system as a wave relayed by the cells. The cells are chemotactically competent and move against the traveling wave of cAMP towards the aggregating centers. In the centers the cells form a mound, where they further develop and differentiate into a slug and eventually into a fruiting body who releases spores which can become myxamoebas; thus completing the life cycle of the amoeba.

The chemical signaling part of the process presents spirals and target centers, these structures are characteristic of both oscillatory and excitable systems, and are present in many systems in nature such as the cardiac muscle [2], calcium waves in oocytes [3], NADH waves in glycolysis [4], among others. Other important system where such structures are observed is the chemical reaction known as Belousov-Zhabotinsky [5, 6], which has become a model experimental system to study spirals and target waves.

A distributed oscillatory system would produce bulk oscillations if it is well stirred, that is, the whole system oscillates synchronously. If it is not stirred and presents some inhomogeneities, such as scratches on the container surface or dust particles, concentric circular waves will appear [7, 8]. We refer to these structures as target patterns. Once a wave breaks, it creates a spiral wave, a persistent structure with topological charge[9, 10].

Many models have been proposed to describe the cAMP waves in D. d.[11, 12, 13], from qualitative excitable models to models derived from the chemical reactions of production of cAMP from ATP. Each model has their own particular way of breaking the homogeneity and producing waves. Here we present the properties of the model proposed by Martiel and Goldbeter [13, 14] when cell inhomogeneity is introduced. We show through numerical simulations that cell inhomogeneity produces naturally stable target patterns which are centered in areas of higher cell density and relayed in the areas of lower density. Therefore, the system behaves both as oscillatory and excitable depending on local density. The produced waves are shown to be of the trigger kind, in contrast to phase waves (or pseudo-waves, without chemical transport), through their characteristic dispersion relation [15]. We also show that wave speed is density dependent. It has been shown in other chemotaxis models that this dependence is enough to produce the wave instability responsible for cell streaming once cell movement is added[16, 17]. Our model reproduces that characteristic streaming pattern of this biological system. Finally we show that a form of degradation which scales with local cell density is needed to have streaming at higher local densities, this degradation can exist in the form of phosphodiesterase bounded to the cell membrane.

Results

Our work combines a cellular automata scheme with the model proposed by Martiel and Goldbeter [13] and extended by Tyson [14] where three fields are described. The extracellular concentration of cAMP is represented as $\gamma$, the intracellular concentration as $\beta$, and $\rho$ represents the percentage of active receptors on the cell membrane. The receptors state changes between an active and an inactive state depending on the cAMP concentration at which they are exposed through the functions $f_{1}$ and $f_{2}$. The intracellular concentration of cAMP $\beta$ increases through production $\Phi(\rho,\gamma)$ and decreases through internal degradation $k_{i}$ and transport towards the extracellular media $k_{t}$. Finally the external concentration of cAMP can diffuse through the system, gets degraded by phosphodiesterase $k_{e}$ and increases due to transport from the internal media $k_{t}$.

We created a grid in our 1-D or 2-D system, dividing it in segments of size $r=10\leavevmode\nobreak\ \mu$m (or respectively squares of area $r^{2}$) and assigned random locations for the cells, each cell had their own value for the intracellular and membrane bound variables $\rho_{i}$ and $\beta_{i}$, while $\gamma(x,y)$ was assigned to a square in space. The used equations look as follows

[TABLE]

[TABLE]

where $s=q\sigma\alpha/(1+\alpha)$ measures the production intensity, $(x_{i},y_{i})$ correspond to the Cartesian coordinates of the $i$-th cell, and $H$ is an indexing function such that $H(i,x,y)=1$ if $x_{i}\in(x-r/2,x+r/2)$ and $y_{i}\in(y-r/2,y+r/2)$ and [math] otherwise. In this way $\gamma$ will increase on a grid space if there is a cell producing cAMP in that location and only diffuse and be degraded if there are no cells in that space. All used parameter are selected following Lauzeral et al. [18] for their good agreement with experimental measures and are indicated in Table 1. We keep $k_{e}$ the external media degradation rate as our control parameter.

The original homogeneous system studied by Tyson et al. [14] presented different regimes as $k_{e}$ was increased. In our set of parameters for low $k_{e}$ the system has one steady state which is stable. At $k_{e}^{*}\approx 4.3$ min-1 the system undergoes a Hopf bifurcation and a stable limit cycle appears. This is the oscillatory regime of the system. Upon further increasing $k_{e}$ two new steady states appear through a saddle-node bifurcation at $k_{e}^{\dagger}\approx 7.74$ min-1, an unstable one and a stable but excitable one. At approximate the same time the limit cycle destabilizes, resulting in the excitable regime of the system. We set our parameter $k_{e}$ so that the system is in the oscillatory regime. For a detailed description of the different regimes present in this system refer to our previous work [19].

Oscillatory clusters

Numerical simulations of a small cluster of very closely located cells producing cAMP surrounded by buffer media without cells, reach a stable steady (non oscillatory) state due to cAMP diffusion to its surroundings. This state is shown in Fig. 1 a). To approximate this solution we calculate separately the area with cells and the area without them. In the area without cells the system reduces to

[TABLE]

which in 1-D has a time independent decaying tail as solution. If the cluster size is $2L$, with the cells located in $x\in(-L,L)$ the decaying tail takes the form

[TABLE]

where $A$ is chosen to fulfill the boundary condition at $x=\pm L$, which are continuity, $\gamma(\pm L^{-})$=$\gamma(\pm L^{+})$, and continuity of the derivative, $\partial_{x}\gamma(\pm L^{-})$=$\partial_{x}\gamma(\pm L^{+})$.

For the cAMP concentration inside the cell cluster we used the fact that simulation results showed that the $\gamma$ values inside the cluster are small compared to the values of the cAMP waves, which justifies making an approximation of the production function $\Phi$ for small values of $\gamma$, we therefore take the polynomial approximation $k_{t}\beta/h\approx a_{0}+a_{1}\gamma+a_{2}\gamma^{2}$ and the system reduces inside the cluster to

[TABLE]

It is worth mentioning that the polynomial approximation values were not taken from the Taylor expansion of the production function for $\gamma\ll 1$, but rather from a least square fitting of the production function. The fitted values for our set of parameters (independent of $k_{e}$) are listed in Table 1. Equation (2) just by itself is invariant to space translation, even though this is not true for the entire system, this invariance produces a conserved quantity, which we will call energy and can be used to calculate the shape of the cAMP accumulated around the cluster. We can rearrange the equation to show its energy conservation more clearly

[TABLE]

and the total energy $E$ can be found by integrating both sides by $\gamma$. Therefore $E=V+(\partial_{x}\gamma)^{2}/2$ is conserved and fixed by the boundary condition. To calculate the shape of the solution we need the energy value, but since the decaying tail has a free parameter we need to introduce a second free parameter and then match the two solutions. Therefore, we leave the maximum value of $\gamma$, $\gamma_{M}$ as a free parameter and we numerically calculated the solution using

[TABLE]

For every $\gamma_{M}$ there is a size $L$ such that it fits the boundary condition to match the decaying tail

[TABLE]

an example of the relation $\gamma_{M}$ vs $L$ is shown in Fig. 1 b). From there, it can be seen that it exists a maximum length for which the static solution exists. For those values of $L$ where two possible values of $\gamma_{M}$ exist, the system chooses the one with smaller $\gamma_{M}$. Indeed linear stability analysis showed that the smaller solution is stable, while the bigger one is unstable. Following this procedure we arrive at two results. First the cluster size values at which a stable, time independent solution exists, and second an approximation of this solution when it exists, see for comparison Fig. 1 a) where both the numerical solution and theoretical prediction are plotted with excellent agreement.

The predicted maximum cluster size for different values is compared to the biggest clusters found through numerical simulations in one and two dimensions in Fig. 1 c) with good agreement for small $k_{e}$. We believe the difference in agreement for large $k_{e}$ comes from the potential $V$ approximation, which is valid only for low values of $\gamma$, since for larger $k_{e}$ the values of $\gamma$ are higher, the approximation is not so good, but still manages to catch the general behavior of the system. For 2-D simulations, cells where located in $r<L$ with $r^{2}=x^{2}+y^{2}$. We see that bigger clusters can be maintained in a low steady state when the degradation is increased, this is expected since for fixed $L$, $\gamma_{M}$ diminishes with increasing $k_{e}$.

If the cell group is bigger than its critical size it oscillates at a frequency that is size dependent, see for example Fig. 2 a), also depicted in Supplementary Video S1, where one oscillating cluster is shown. As the cell group becomes larger, the frequency increases converging towards the limit cycle frequency. We measured this value for different cluster sizes once the cluster frequency converged to measuring precision. This relation is shown in Fig. 2 d).

At small sizes the cluster oscillation is synchronized (bulk), but as it gets bigger ($L\gtrapprox 0.19$ mm) a wave develops. This wave starts at the center of the cluster and moves towards the sides increasing its amplitude as it travels, as can be seen in Fig. 2 b) and c). Compare to Fig. 2 a) which shows a smaller cluster which oscillates synchronously. This effect can also be seen in Supplementary Video S2 which in turn can be compared to Supplementary Video S1.

Target patterns in random distributions of cells

When we assigned random locations to the cells, the system showed different regimes as we increased $k_{e}$. For low values of $k_{e}$ the system oscillated mostly synchronously, see Supplementary Video S3 and Fig. 3 a) for a space-time plot. For high values of $k_{e}$ the system was incapable of oscillating by itself and reached a steady state of low cAMP (see Fig. 3 c)).

A more interesting behaviour was observed for intermediate values of $k_{e}$. There, the areas of higher local density acted as oscillatory centers, while the rest of the system relayed the emitted waves. This means that the lower density areas although not capable of oscillating by themselves are capable of producing enough cAMP to maintain the wave and avoid its complete degradation.

This heterogeneity in the system response, seemingly both oscillatory and excitable, appeared even though the same parameters were used trough the whole system, the heterogeneity being given by the cell distribution.

The appeared target centers have a range of frequencies and are stable. They also interact with each other, thus higher frequency centers dominate over the lower frequency ones. These characteristics match the observation in D.d. patterns [20]. A typical example of our simulation results can be seen in Fig. 3 b) and in Supplementary Video S4.

The range of degradation values $k_{e}$ in which centers can be observed depends on cell density as is shown in Fig. 4a). Under the shaded area bulk oscillations like in Fig. 3 a) are observed, and above the shaded area spontaneous center do not appear, like in Fig. 3 c). At higher cell densities, higher degradation rates are required to observe spontaneous centers, which is consistent with the idea that phosphodiesterase is produced and released to the external media by the cells, therefore a higher concentration of cells should produce a higher concentration of phosphodiesterase and therefore a faster degradation. Nevertheless, the amount of maximum degradation does not scale linearly with density, but the size of the existence range decreases with higher density.

The transitions between synchronized oscillations, center creation, and low steady state are not bifurcations but rather smooth and of a stochastic nature. That means, for example, that for high values of $k_{e}$ a small number of centers may still appear in some simulations if by chance big clusters of cells appear together. The shaded area of Fig. 4 a) indicates where centers spontaneously appear in more than $75\%$ of the cases. See Methods for a description of the boundary calculations for Fig. 4 a).

To understand the nature of the traveling waves we calculated their dispersion relation, i.e. the shape of the velocity-period curve. We used the homogeneous distribution simulations mentioned in Methods to study more controlled waves. The waves are generated perturbing the system with a period $T$ and the velocity is calculated following each individual peak once they passed the initial transient. These dispersion relations are presented in Fig. 4 b). It can be seen that they present the characteristic shape of trigger waves [15], therefore they are actual waves that produce transport of chemicals in contrast to pseudo- or phase waves which do not involve chemical transport and have a smaller dependence on diffusion. Here, diffusion is necessary to activate the next cell and with that relay the wave, another characteristic of trigger wave behaviour.

If the homogeneous system is perturbed with higher frequencies than those depicted in Fig. 4 b) it will not react with a 1:1 response, but instead it will propagate the waves with a lower frequency, usually located in the elbow (strong curvature) area of the curve.

Furthermore, we measured the frequency of the signaling centers in 1-D simulations for a fixed cell density of $4\cdot 10^{5}$ cells/cm2 (0.4 mono-layer) as a function of degradation rate. We measured lower center frequencies at higher degradation rates, as shown in Fig. 4 c).

Streaming

One common test for the waves generated in a D.d. model is if the system is capable of showing streaming once cell movement is added. Streaming is a characteristic feature of the intermediate state of D.d. on their way towards the mounds. It consist on the alignment of the cells in a head to tail manner, thus displaying long ramified lines of higher density that spread radially from the aggregating center.

To test this feature, we added cell movement in a simple toy-like manner. The cells would move at a constant speed $v_{c}$ if they sense a cAMP gradient bigger than a minimum $\partial_{x}\gamma>\Delta\gamma_{c}$ and if their percentage of active receptors is bigger than a cutoff $\rho>\rho_{c}$. The first rule avoids movement due to random noise and the second avoids the back-of-the-wave problem [21, 22]. The back-of-the-wave paradox consists in taking into account that D.d. only moves in the first half of the wave, ignoring the gradient of the decreasing cAMP in the second half which would move the cell in the opposite direction. By adding a minimum $\rho$ to allow the cells to move, they effectively move only in the first part due to the desensitization produced by the passing wave. The cells continue moving as long as these both conditions are met. Initially we did not allowed cell superposition, therefore a cell would only move towards a different space if the arriving location were empty. The updated cell positions were calculated each time step using a forward Euler scheme, after the $\gamma$, $\beta$, and $\rho$ fields were calculated. All used parameter are listed in Table 1.

We performed numerical simulations in 2-D where we allowed cell movement after a pattern was already established ($t\approx 50$ min). Images of a typical simulation are presented in Fig. 5 where the streams can be clearly observed, thus recovering the expected behavior. The full simulation can be seen in Supplementary Video S5. It is also worth mentioning that after cells start to move some wave fronts may break and produce spirals.

Cell superposition and bounded phosphodiesterase

As a way of perfecting our simple toy model for movement we included the possibility of cells superimposing in the same location, with a maximum of 5 cells per grid point. In order for the system to continue presenting centers and continuous cell streaming, additional degradation is needed. That is, with the previously used constant degradation the centers and aggregation streams start to break once cell clumps appear due to movement. See Fig. 6 in comparison to Fig. 5 where target centers and stream lines quickly break apart. A simulation showing this behaviour is presented in Supplementary Video S6.

Many experimental studies point to the existence of phosphodiesterase bounded to the cell membrane [23, 24]. Adding this extra bounded phosphodiesterase $k_{eB}$ solves these problems and the system behaves as before, producing centers which persist after movement (see Supplementary Video S7). Therefore, we conducted simulations where there is a constant degradation present on the system $k_{eU}$ which represents the unbounded, free phosphodiesterase; and a bounded degradation $k_{eB}$ which exists attached to the outside of the cell membrane and therefore exists only on the occupied spaces and in an amount dependent on the number of cells in that location. The new equations are then

[TABLE]

The results obtained with these equations are similar to those shown in Fig. 3, that is, we distinguish three types of behaviour depending on the degradation rates. The degradation combinations for which the system still presents spontaneous centers are shown in Fig. 7. Unlike in the case with no bounded phosphodiesterase, the variance of the normal distribution decreased with the bounded degradation rate, thus making the transition from centers to low steady state much sharper at higher $k_{eB}$.

We see that there is no fixed ratio of bounded/unbounded phosphodiesterase such that the system produces centers, but rather there is a whole range of values, depending on the system density, thus allowing the cells to have some variability and still present the same behaviour. It is noteworthy that as the bounded degradation increases ($k_{eB}$), $k_{eU}$ decreases, as does the range of possible unbounded values $k_{eU}$, and the upper transition becomes sharper, since the range of values where $50-75\%$ of simulations presented centers diminishes (length of the upper errorbars). Thus the cells can have less variability in $k_{eU}$ as $k_{eB}$ increases.

Discussion

Over the past years several models have been used to describe the patterns observed in D.d. Excitable models are capable of sustaining spiral and trigger waves, but in order to generate a target pattern they require either an oscillatory center or spontaneous/random firing from the cells. On the other hand, oscillatory models require some sort of inhomogeneity to be added to the system to produce stable centers and avoid bulk oscillations. It has been shown that increasing the cell parameters along a developmental path or adding random firing decides the location of the observed patterns in an artificial way [25, 26, 27, 28] which contrast with experimental observations showing that oscillation is a collective effect instead of the work of some specialized cells, where even the cells composing the oscillating center move continuously in and out of the signaling center [29, 30]. We suggest that a more simple mechanism is also in play to produce centers for densities less than a mono-layer: the inhomogeneous cell distribution in the system is enough to create stable emitting centers. Our model matches the observation of Lee et al. [31] who showed that at lower densities (below mono-layer density $\approx 10^{6}$ cells/cm2) target centers dominate the observed patterns, while high above confluency, spirals appear. Our simulations show that below confluency and for a fixed value of $k_{e}$, the number of target centers increases with cell density until the transition to bulk oscillations occur (see Fig. 4a)). A recent experiment by Ohta et al [30] showed that at densities of the order of $1.25\cdot 10^{5}$ cells/cm2 (0.125 mono-layer) centers can be generated by groups of roughly 13 synchronized cells in an area of $100\times 100\leavevmode\nobreak\ \mu$m2. These results match our observation with $k_{e}=1.56$ min-1 where a group of approximately 18 cells produce an oscillatory center, see Fig. 8 and Supplementary Video S8.

At higher cell densities when the pattern evolves towards spiral creation, we expect this mechanism to be less relevant, although our model did show occasionally spirals when the wave front broke due to high inhomogeneities in cell density (see Supplementary Video S9). Above mono-layer densities it has been shown that the cell development inhomogeneity [18, 32] and the strength of the feedback loop on excitability [33] play a prominent role in the spiral creation process.

We stress the simplicity of this mechanism which does not make big assumptions about the inside state of the cells, but simply requires an external degradation mechanism to be present in the media. Small groups of isolated cells reach a steady state of low cAMP concentration. We showed that this stable solution exist for small clusters, where the system has two solutions, a stable and an unstable one. Once the cluster size increases over a critical value, the cluster is no longer stable and starts to oscillate with a size dependent frequency, therefore bigger clusters have a higher frequency which allows them to dominate over smaller clusters. In simulations with randomized cell locations the centers are not necessarily large clusters of consecutive cells, but rather just areas of a higher local density, where some clusters are usually separated by small distances and act together as a big cluster. Even though this does not exactly match the scenario described in the Oscillatory Clusters section, it does provide an insight into why these centers spontaneously appear. Our model is therefore, a collective effect where groups of cells behave oscillatory and individual cells are excitable to external cAMP pulses, consistent with experimental observations [29], furthermore the lower cAMP level at lower densities is fundamental for their excitable behavior, displaying an oscillatory behaviour if the cAMP levels are artificially risen[29].

Having a constant degradation across the entire system is an oversimplification of the real setup where the cells produce different amounts of phosphodiesterase [34] and it even changes over time [35], but since the centers exist over a range of degradations (see Fig. 4a)) we believe that the observed effect is robust to individual differences of phosphodiesterase production in cells, an effect which is also dampened by phosphodiesterase diffusion. Further modeling including phosphodiesterase diffusion and production is necessary to explore these possibilities.

The target patterns showed to have all the main features of those observed in experiments, like producing trigger waves, having different frequencies, and higher frequency centers taking over lower frequency ones. The system showed a range of degradation values at which it presents centers, this range decreased with increasing density, while the maximum value scaled non linearly with density. For a fixed cell density increasing the degradation rate decreased the average center frequency (see Fig. 4c)), consistent with measurements in flow chamber experiments [29].

Once cell movement was introduced, the system showed streaming when aggregating. We attribute this behaviour to the velocity dependence on density (see Fig. 4 b)), as it has been shown before in similar descriptions [16, 17, 36]. In those systems, as in ours, the cAMP wave travels faster on the areas with higher cell density thus producing a wave which is not perfectly circular. This shape leads to an aggregation which is not radially uniform, thus producing streams. In this model the velocity dependence on density comes from the wave speed difference between spaces with cells, where new cAMP is produced, and empty spaces where speed is given only by diffusion, thus on average higher density areas relay waves faster than lower density ones, producing streaming. This simple aggregation model shows that intracellular localization of cAMP or ACA is not necessary to explain streams, along with other mechanisms such as cell-cell adhesion, chemotactical memory [37], and directional rectification [38], that, although are present in the experimental system, do not seem to be fundamentally necessary for stream formation.

The final section highlights the role of a cell bounded form of degradation. This comes into relevance once the cells have started to move and form small clumps. In order for these clumps to not disrupt the wave propagation process, a degradation that scales with density is necessary. The existence of this type of phosphodiesterase is still a matter of discussion in the experimental community[39, 34, 24, 40] and even though we believe our model provides arguments in favor of the existence of both types, the effect of locally increasing degradation can also be explained by a local accumulation of phosphodiesterase around the cells possibly due to the recent release of PDE to the extracellular media that has not yet diffused sufficiently.

In conclusion, we have presented a scenario of creation of target patterns in a model describing pattern formation in D.d.. This scenario does not require specific stages along the cell developmental path, nor requires the introduction of random firing, both of which preselect the location of target centers. Instead, in this description wave centers appear spontaneously in areas of higher cell density. We have shown that introducing cell inhomogeneity in the Martiel-Goldbeter model creates naturally target patterns that are stable and capable of producing waves that fill the whole system. We characterized these centers, which correspond to areas of higher cell density. The size of the minimum cluster required to produce a center increases with the degradation rate of the system. If the cluster size is not big enough to sustain autonomous oscillations, it reaches a low steady state whose amplitude decreases with the degradation rate. These smaller clusters, or areas of lower local density can nevertheless be excited, thus allowing waves to propagate through the system. Therefore, this scenario reproduces the large scale organization displayed by D.d. populations.

Methods

The system was simulated using a fourth-order Runge-Kutta scheme with Merson correction [41] and an adaptive time step. We used finite differences and a 5 points laplacian in 2-D and 3 points in 1-D. For 1-D simulations $dx$ was selected equal to the cell size $r=10\leavevmode\nobreak\ \mu$m. To speed up calculations for most 2-D simulations, $dx=5r$ was selected, with the results being confirmed by smaller grid sizes.

To calculate the upper boundary of Fig. 4 a) and Fig. 7 we simulated the system a hundred times for each set of parameters and recorded how many simulations presented centers. We then fitted a gaussian distribution around the boundary and from this fitting extracted the expected value for a $75\%$ of center appearance (marked as a circle) and for a $50\%$, marked as the end of the errorbar. In Fig. 4, where only unbounded phosphiesterase was present, the fitted gaussian distributions presented the same variance, differing only in their mean. In Fig. 7 the variance lowered with increasing $k_{eB}$.

To calculate the lower boundary we used a extreme case consisting in homogeneously distributed simulations, where cells had a constant distance between them. This is the most spread out possible distribution for the cells, thus avoiding clusters as much as possible. For example, if the density is $5\cdot 10^{5}$ cells/cm2 (0.5 mono-layer) the distribution would be one cell, then 0.01 mm of empty space, then another cell and so on. In this homogeneous system if the degradation is too small the cells show bulk oscillations. As we increased the amount of phosphodiesterase the bulk oscillations ceased and the system reached a stable steady state of low cAMP concentration. The lower boundary of Fig. 4 a) and Fig. 7 indicates the maximum amount of phosphodiesterase such that the homogeneous system shows bulk oscillations, therefore at that degradation rate or higher any inhomogeneous cell distribution will produce target centers.

Data availability statement All data used in this study will be made available upon request.

Acknowledgements

The authors thank E. Bodenschatz, A. Bae, and V. Zykov for fruitful discussions. E.V.H. thanks the Deutsche Akademische Austauschdienst (DAAD), Research Grants—Doctoral Programs in Germany. A.G. acknowledges MaxSynBio Consortium, which is jointly funded by the Federal Ministry of Education and Research of Germany and the Max Planck Society.

Author contributions statement

E.V.H. and A.G designed the research, E.V.H. performed the numerical simulations and theoretical analysis, E.V.H. and A.G wrote and reviewed the manuscript.

Additional information

Competing interests The authors declare no competing interests.