Segmentation-clock synchronization in circular-lattice networks of embryonic presomitic-mesoderm cells

1.   Introduction

Synchronization is an emergent phenomenon prevalent in many oscillatory living systems, and it is crucial in numerous processes, either for maintaining life or to cause diseases ^[1,2,3]. One of these processes is somitogenesis. Let us briefly review the main embryonic stages leading to somitogenesis (the formation of somites). The gastrulation stage is characterized by diverse cell migration and differentiation steps that give rise to three different layers: ectoderm, mesoderm, and endoderm. From these three layers, all embryonic tissues will derive ^[4,5]. Subsequently, the mesoderm is divided into the lateral mesoderm, the intermediate mesoderm and the presomitic mesoderm (PSM). The latter is then segmented during somitogenesis, giving rise to blocks of cells called somites ^[6,7]. Somites are the precursors of various tissues of the back as: ribs, vertebrae, skeletal muscle and dermis ^[5,8]. Disruption of somite formation can result in the occurrence of syndromes and disorders that include abnormal vertebral segmentation, such as spondylocostal dysostoses ^[9]. Somites are progressively formed in pairs from the anterior end of two rods of PSM tissue, which lie on both sides of the caudal neural tube, and join at the posterior end of the embryo (a region known as the tail bud) ^[10]. A characteristic feature of somitogenesis is its periodicity, which is controlled by a complex gene network termed the segmentation clock ^{[11,12,13,14]}.

The segmentation clock drives the periodic expression of a number of genes in PSM cells, with a periodicity that matches somite formation and is species specific: $\sim\, 90$ min in chick, $\sim\, 120$ min in mouse and $\sim\, 30$ min in zebrafish ^[12,15,16]. Although the detailed clockwork of the segmentation clock is quite complex, the basic mechanism is that of a self-inhibiting gene. This negative feedback loop, along with the delays involved in transcription, splicing, translation, post-translation modifications, etc., are responsible for giving rise to periodic oscillations in the concentration of the encoded proteins ^[17]. Figure 1 summarizes this basic gene network for the zebrafish segmentation clock.

Despite the many recent advances unveiling the segmentation clock clockwork in various species, it was not clear until very recently whether individual PSM cells are able to maintain a rhythm of gene expression on their own, or they rely on the interaction with surrounding cells to achieve it. Webb et al. ^[18] removed individual PSM cells from zebrafish embryos and recorded, by means of a fluorescent reporter, the expression level of gene her1 (a known segmentation-clock gene) over time. Webb et al. were able to show that individual cells are capable of maintaining a rhythmic pattern of her1 expression on their own. However, they do it in an intermittent way, and when comparing the her1 activity patterns in individual cells with the patterns of larger segments of zebrafish tissue, Webb et al. found that rhythms in individual cells are slower and have a less precise timing than those in the tissue. We speculate that, although individual PSM cells are capable of showing autonomous oscillations, the precision and timing needed for somitogenesis to occur correctly is an emergent phenomenon of groups of PSM cells that depends on their ability to synchronize.

Single PSM cells interact with adjacent cells through the Delta-Notch signaling pathway, and this interaction allows them to synchronize ^[19,20,21]. The Delta-Notch signaling pathway in somitogenesis works as follows. Notch is a transmembrane protein composed of three segments: the Notch extracellular domain (NECD), the transmembrane domain (TM) and the Notch intracellular domain (NICD) ^[22,23]. NECD functions as a receptor to various ligands, but in PSM cells it binds to ligands of the Delta family. Delta is a transmembrane protein whose expression is controlled by the segmentation-clock regulatory mechanisms; thus, its concentration oscillates periodically ^[21]. Once the Delta ligand from one cell binds to the Notch receptor of a neighbor cell, the NICD is released and is capable of reaching the cell nucleus to act as a transcription factor (activating transcription) for target genes that form part of the segmentation clock clockwork ^[24,25]. A representative diagram of the mechanism described above is shown in Figure 2.

Numerous works have studied synchronization in PSM cells. For instance, Lewis ^[17] showed by means of a mathematical model that two adjacent cells interacting via the Delta/Notch pathway can synchronize, and studied the effect of time delays upon synchronization. Tiedemann et al. ^[26] studied the influence of time delays on the synchronization of a network of PSM cells located in a regular square lattice. Riedel-Kruse et al. ^[27] analyzed the clock's synchrony dynamics by varying strength and timing of Notch coupling in zebrafish embryos with techniques for quantitative perturbation of gene function, and developed a physical theory, based on coupled phase oscillators, which explained the observed onset and rescue of segmentation defects, the clock's robustness against developmental noise, and a critical point beyond which synchrony decays; Horikawa et al. ^[28] examined the response of a normally oscillating clock in zebrafish to experimental stimuli using in vivo mosaic experiments and mathematical simulation, and found that intercellular coupling has a crucial role in minimizing the effects of noise to maintain coherent oscillation; and Uriu et al. ^[29] demonstrated that random cell movement promotes synchronization of the segmentation clock. Despite all of these advances, some interesting questions like: how the size of a network of PSM cells, the number of interactions per cell, the strength of cell-cell interactions, and the inherent variability in the PSM cells, affect segmentation-clock synchronization, remain unanswered. The present work is aimed to tackling these questions from a mathematical modeling perspective.

2.   Methods

2.1.   Mathematical model

The model here employed to study PSM-cell synchronization is based on the zebrafish segmentation-clock clockwork—see Figure 1. In this system, the commonly accepted explanation for the origin of single PSM cell oscillations is the auto-repression of genes her1 and her7, after the corresponding proteins form dimers with Hes6—another member of the basic helix loop helix (bHLH) protein family that plays the role of a transcriptional repressor and has no oscillatory expression ^{[29,30,31,32]}. These three proteins assemble in dimers in any combination, but only homodimers Her1/Her1 and heterodimers Hes6/Her7 have a repressor effect on the expression of her1 and her7. In ^[33], Schroter et al. employed a mathematical model to study the dynamics of the Her1, Hes6 and Her7 proteins at the single cell level. The interactions they considered are shown in Figure 1. By taking into account these interactions, making some plausible biochemical assumptions, and normalizing the state variables, Schroter et al. developed the following system of delay differential equations (DDE) to model the behavior of proteins Her1, Hes6, and Her7 within a cell:

In the above equations, $h_i$ ($i = 1, 6, 7$) stand for the normalized protein concentrations, $\gamma_i$ denote the corresponding maximum transcription rates, $\delta$ is the effective dimer degradation rate, and $\tau_i$ correspond to the necessary time to yield a functional protein from transcription initiation to post-translational modifications. Details on how these equations are obtained can be found in ^[33]. Notice in the above equations that the terms accounting for the degradation of $h_1$, $h_6$ and $h_7$ are equal to the normalized concentrations of the corresponding chemical species. This is because Schröter er al. ^[33] assumed that all proteins degrade with the same rate, and normalized time to the inverse of the protein degradation-rate constant.

Since we are interested in studying synchronization in networks of PSM cells, we have to consider cell-to-cell interactions. Hence, we took the model by Schröter et al. as a starting point, and expanded it by taking into account that coupling between PSM cells occurs via the Delta-Notch signaling pathway. The gene coding for protein Delta is repressed by the same dimers that repress genes her1 and her7 ^[34]—see Figure 2. Therefore, after making the same simplifications as in ^[33], a differential equation for the dynamics of protein Delta can be written as

where D denotes the concentration of Delta protein, while $\gamma_D$ and $\tau_D$ are defined as in equations (2.1)-(2.3). Observe that the first term in the right hand side of Eq (2.4) has the same functional form as the corresponding terms in Eqs (2.1) and (2.3). That is, we have assumed that the gene coding for protein Delta is regulated by homodimers Her1/Her1 and heterodimers Hes6/Her7 in the same way as genes her1 and her7. The second term in the right hand side of Eq (2.4) corresponds to the assumption that protein Delta degradation is driven by exponential decay. We further assumed that protein Delta degrades with the same rate as proteins $h_1$, $h_6$ and $h_7$. From this assumption and the time rescaling introduced by Schröter er al. ^[33], the term accounting for the degradation of protein Delta is proportional to the corresponding protein concentration.

Equations (2.1)-(2.4) constitute a model for the dynamics of the zebrafish segmentation-clock, including protein Delta, at the single cell level. To account for a network of PSM cells that interact via the Delta-Notch pathway, this model can be amended as follows:

In the above equations, subscript $i$ refers to the $i$-th cell in the network, while the sum in the argument of function $g(\sum_jD_j)$ is performed over all $j$ cells that interact with the $i$-th cell via the Delta-Notch pathway. Moreover, function $g(\sum_jD_j)$ stands for the effect that Delta proteins in the interacting cells have upon the transcription rate of the $i$-th cell genes. Since Delta proteins in neighboring cells promote the release of Notch intracellular domains (NICDs) within the $i$-th cell, while NICDs act as activators of genes her1 and her7, we assume that function $g(\sum_jD_j)$ can be modelled via the following increasing Hill function ^[35,36]:

where $k_D$ is the corresponding half saturation constant.

The way Delta regulation has been introduced in Eqs (2.8) and (2.9) is consistent with the assumption that protein Delta enhances the expression of genes her1 and her7 by binding to a different DNA segment than dimers Her1-Her1 and Hes6-Her7, and that they act independently. Although other regulatory functions can be pictured, we consider that this is a plausible one.

Despite some evidence that Notch levels are also dynamic in the segmentation clock, we have assumed, for the sake of simplicity, that Notch levels remain constant and that only Delta levels change in time.

2.2.   Parameter estimation

The model parameter values considered in this work were taken from ^[33,34]. All parameter values are tabulated in Table 1, after rescaling them to account for the normalization process.

Parameter $k_D$ is a control parameters whose value is changed to investigate its influence on the system behavior.

2.3.   Order parameter

To test whether two or more oscillators are synchronized, we employed the order parameter $R$ defined in ^[37]. This parameter, which quantifies the distribution of phases of a group of oscillators, is defined as follows:

where $Var (M)$ refers to the variance of the mean time-series, averaged over all oscillators; while $\overline {Var (b_i)}$ is the average of the individual time-series variances. In the definition given by Eq (2.10), the numerator measures the distribution of oscillator phases, while the denominator is a normalization factor that makes the order parameter range between 0 (no synchronization) and 1 (perfect synchronization, with all oscillators in phase) ^[38].

2.4.   Numerical simulations

The model equations were numerically solved in $\texttt{Python}$ by means of $\texttt{SciPy}$'s $\texttt{PyDelay}$ library, which uses the Bockagi-Shampine method for solving delay differential equations.

Notice from Table 1 that $\tau_7 = 1$. This means that, in the non-dimensional model, time has been normalized to the delay associated to transcription, mRNA splicing, translation, and post-translational modifications of protein Her7, which is about 18 minutes ^[17]. Taking this into consideration, all simulations were run for 5500 normalized units of time, but only the time series corresponding to the last 500 units of time were taken to test synchronization. In this way, we make sure that the system has reached a stationary behavior, independent of the initial conditions.

In our simulations, we took randomly generated initial conditions and, in some cases, randomly generated parameter values. Regarding the initial conditions, all the variables of each network cell were set to an initial value randomly selected from a uniform distribution in the interval $[0, 1]$. That is, all variables corresponding to a given cell have the same initial value, while all cells have different initial states. In this way, we ensure that all network cells have different initial phases. Since the outcome of each simulation is different, due to the various sources of randomness, we performed 10 independent simulations for every network condition to be able to statistically analyze the obtained results.

3.   Results

We are interested in understanding how parameters like: network size, connectivity level, interaction strength, and cell-to-cell variability, affect synchronization in a network of PSM cells. To tackle these questions, we considered very simple and idealized networks, which are known as circular lattices—see Figure 3. The vertices in these networks represent cells which are located on a ring. Furthermore, each vertex is connected by edges (that represent mutual interactions) to $m$ vertices ahead and $m$ vertices behind. We say that the connectivity level of a network like this is $m$. For the sake of simplicity, we only considered networks with and odd number of vertices. Thus, the maximum possible connectivity level of a network with $n$ vertices is $(n - 1)/2$.

Although circular lattices are unrealistic representations of real PSM cell networks, both of them share in common that each vertex (cell) is connected with a reduced number of neighboring vertices (cells). Recall that interactions among PSM cells only occur when there is physical contact between them. In our opinion, this resemblance makes circular lattices suitable to study how network size, connectivity level, interaction strength, and cell-to-cell variability, affect the synchronization of PSM cell networks, without having to deal with other phenomena, like border effects.

3.1.   Effect of network size and cell connectivity on synchronization

To investigate how the connectivity level of a network of PSM cells affect synchronization, we simulated the time evolution of networks consisting of 21 vertices (in which all vertices represent identical PSM cells) and different connectivity levels, with parameter values as given in Table 1 (these parameter values allow individual cells to spontaneously oscillate), while the parameter corresponding to the strength of interaction between cells, $k_D$, was set to $k_D = 1$. In each simulation, the initial conditions were set as described in the Numeric Simulations subsection. In Figure 4 we present the results of sample simulations with connectivity levels equal to 1, 3, and 6. Observe that, as the connectivity level increases, clusters of synchronized cells start to appear, and that all cells synchronize at connectivity level equal to 6. The phenomenon in which certain clusters inside the network show synchronization is known as cluster synchronization, and is a current active field of study; see for instance ^{[39,40,41,42,43]}. However, to the best of our knowledge, this is the first time it is reported in the segmentation clock literature.

To quantify and expand the results in Figure 4, we considered all possible connectivity levels for such network, performed 10 simulations for each network condition, and computed the average order parameter as a function of the connectivity level. The same was done for networks of sizes 11 and 41, to study the influence of network size on synchronization. The obtained results are summarized in Figure 5. Observe that increasing the network connectivity level (recall that in a network with connectivity level $m$, each cell interacts with $2m$ other cells) facilitates synchronization of all cells. This behavior does not seem to depend on the network size. Interestingly, synchronization of all the cells in the network is always achieved when the connectivity level is equal to or larger than 6. That is, when each cell in the network interacts with at least its 12 nearest neighbors (6 ahead and 6 behind along the circular lattice). This further suggests that there is a minimal network size that allows synchronization of all its cells. For the present parameter set, the network has to consist of at least 13 cells to achieve synchronization.

As mentioned before, we decided to investigate the synchronization of segmentation clocks in symmetric circular-lattice networks, due to their simplicity and their having certain properties in common with real networks. As for their simplicity, they allow us to study the effect of network size, without worrying about border effects (like in square lattices); as well as to study the effect the connectivity level of has on synchronization, without dealing with network variability, as in random networks. To test the generality of our results, we repeated the simulations in Figure 5, but considering random networks with the same connectivity level, in which the vertices each vertex is connected to were randomly selected. In the Supplementary Data section we provide a link to a \texttt{Jupyter} notebook containing the results of such simulations. We can see there that the obtained results are quite similar to those from circular lattices, thus proving that they are not specific for circular lattices.

3.2.   Robustness of synchronization to variability of the coupling-strength

We are also interested in studying how the strength of inter-cellular interactions affect synchronization. In our model, inter-cellular interactions are accounted for by function $g(D)$, which stands for the effect that Delta-Notch-mediated interactions with neighboring cells have on the expression level of the segmentation-clock genes. We can see in the definition of $g(D)$—Eq (2.9)—that smaller $k_{D}$ values mean that lower concentrations of Delta proteins in neighbor cells are required to activate the segmentation clock genes. From this, we assume that the strength of inter-cellular interactions is a decreasing function of $k_{D}$.

To investigate how weakening or strengthening cell-to-cell interactions affects synchronization, we performed simulations with a network of size 21 and connectivity level equal to 2, and varied parameter $k_D$ in the range $[0,164]$. This wide range of parameter values allowed us to observe the repertoire of behaviors illustrated in Figure 6. That is, cells are not synchronized when $k_D$ is zero or very close to zero (setting $k_D = 0$ is equivalent to assuming that all network cells are independent because $g(D) = 1$ in such case). As the value of $k_D$ increases, clusters of synchronized cells start to appear, and they increase in size until, eventually, all cells synchronize. Full synchronization holds for a while but, when $k_D$ is so large that cell-to-cell interaction is weak enough, the network cells oscillate with different phases once more. If the value of $k_D$ is further increased, a regime in which all cells oscillate in synchrony but with a negligible amplitude is eventually reached.

We further performed several more independent simulations for the same type of networks. We carried out 10 simulations for each $k_D$ value and computed the corresponding order-parameter mean value. The observed dependence of $R$ on $k_D$ can be appreciated in Figure 7. This figure gives us a clearer picture of the $k_D$ ranges in which the different behaviors illustrated in Figure 6 are observed.

The results in Figures 6 and 7 correspond to a very specific network. To test their generality, we investigated the behavior of networks of various sizes and connectivity levels. The results are summarized in Figure 8. We can observe there that the same repertoire of behaviors is observed, regardless of the network size and connectivity level. For the sake of clarity, we decided not to show the segments where curves $R$ vs. $k_D$ rise to values close to 1 for large $k_D$ values, but this behavior is certainly observed in all cases. Further notice that synchronization can be achieved with larger $k_D$ values as the connectivity level increases. This indicates that increasing the number of interactions per cell compensates the reduction of coupling strength. In other words, more interacting cells are required to maintain a given expression level of protein Delta when the value of $k_D$ is increased.

We are particularly intrigued by the appearance of negligible-amplitude synchronous oscillations in the large $k_D$ regime. This phenomenon, which has been termed "oscillation death" ^[44], is usually present when a group of oscillators have a strong coupling. However, it occurs in our system in the weak coupling regime. We cannot offer an explanation for why this happens, but it may be related to the presence of time delays. In any case, this deserves further investigation.

To summarize the results in Figures 6–8, we elaborated the phase diagram in Figure 9, which shows how synchronization depends on parameter $k_D$ and the network connectivity level. This diagram corresponds to a network of size 21, but qualitatively similar results are obtained for different network sizes. Such a diagram is the result of analyzing the value of parameter $R$ and the amplitude of oscillations when varying $k_D$ in the range [0,164] (in a grid of 18 points equally spaced along the logarithmic $k_D$ axis) for each connectivity level.

We also explored how the above results depend on the various parameter values. We found that they are not noticeably affected, even by quite large parameter variations, except for $\gamma_D$, the maximum transcription rate for protein Delta. In the Supplementary Data section we provide a link to a public repository with additional plots containing the results of these extra simulations. We can see there that increasing $\gamma_D$ favors synchronization. However, increasing the connectivity level has a stronger effect.

3.3.   Robustness of synchronization to parameter variability among PSM cells

After analyzing how the capability of PSM-cell networks to synchronize depends on network size, connectivity level, and interaction strength, we decided to investigate the robustness of synchronization to parameter variability among PSM cells. We performed simulations in which all the parameters of every cell in the network were randomly selected from normal distributions (with mean equal to the corresponding base values reported in Table 1, while $k_D = 1$, and with standard deviations ranging from 5% to 25% of the corresponding base values). Since the largest standard deviation we considered is 0.25 times the mean value, the probability that negative values are sampled is very small. However, to avoid this possibility, we added an instruction in our algorithm such that, whenever a randomly sampled parameter resulted negative, its value was set to zero.

In Figure 10 we show the stationary behavior of individual simulations computed with a network of size 21, connectivity level equal to 4, and different levels of parametric variability. Notice how, in this particular network, larger levels of parametric variability induce cluster, and in one case full, synchronization.

To test the generality of the results in Figure 10, we performed similar simulations, with networks of different sizes and different connectivity levels. We performed 10 independent simulations for each network condition, and computed the corresponding order-parameter mean value and standard deviation. The obtained results are summarized in Figure 11. Observe that parametric variability among cells has a positive effect on synchronization (it increases the value of the order parameter, $R$) for connectivity levels smaller than 6, regardless of the network size. To better illustrate these observations, we present in Figure 12 plots of $R$ vs. the amount of parametric variability, for networks of different sizes and a connectivity level equal to 4. We can also notice in Figure 11 that parameter variability has a deleterious effect on synchronization for connectivity levels larger than 6, since the value of $R$ decreases as the amount of parametric variability increases. These deleterious effects of parameter-variability are attenuated as the network size increases. In summary, the above results suggest that a moderate level of variability among PSM cells may favor synchronization in networks with a connectivity level of 5 or less.

4.   Concluding remarks

We have studied (via a mathematical modeling approach) the synchronization of segmentation clock oscillations in very simple networks of zebrafish PSM cells. Our results can be summarized as follows:

● There exists a connectivity-level threshold below which synchronization is impossible. Such threshold depends on the value of $k_D$. This result implies that very small PSM-cell networks would not be able to synchronize.

● Decreasing the strength of interactions among PSM cells has a deleterious effect upon synchronization. However, this effect can be compensated by increasing the network connectivity level.

● Oscillation death takes place in the extremely-weak-coupling regime.

● A moderated level of variability among network cells can have a positive effect on synchronization, specially at small connectivity levels.

The second result agrees with observations on synchronization studies on networks of other types of oscillators ^[45,46]. Thus, it was somehow expected. Contrarily, as far as we are concerned, the other results are novel and may have interesting implications. Regarding the existence of a connectivity-level threshold for synchronization, one would expect by geometrical considerations that, in a 3-dimensional PSM, each cell would be in contact with about 10 cells, which agrees with the order of magnitude of the threshold we found. In the same line if reasoning, Delaune et al. ^[47] experimentally observed that the average number of interactions per cell is about 8-10. The fact that cell variability enhances synchronization is also interesting, and its biological implications are straightforward (there is no such a thing as two identical cells). We find intriguing that oscillation death takes place in the weak coupling regime, and understanding why this happens would be appealing from the point of view of nonlinear dynamics. We wish to point out that, from a nonlinear dynamics perspective, our work lies in the field of synchronization in time delayed systems, which is a relatively new and unexplored research field; see for instance Refs. ^[48,49].

5.   Supplementary data

Here, we provide links to Jupyter notebooks with additional plots for Sections 3.1 and 3.2.

Additional plots for section 3.1:

https://github.com/JesusPantoja/PSM-Synchronization-Plots/blob/master/PlotSuppVCIR3d.ipynb

Additional plots for section 3.2:

https://github.com/JesusPantoja/PSM-Synchronization-Plots

Conflict of interest

All authors declare no conflicts of interest in this paper.