A mathematical framework for measuring and tuning tempo in developmental gene regulatory networks

Biological systems are intrinsically dynamic. Cells are constantly migrating, differentiating, dividing and dying; collectively, these activities underpin numerous biological processes. This dynamism is especially apparent during embryonic development, during which the precise coordination of cell differentiation plays a pivotal role in the orchestrated formation of distinct tissues constituting the anatomical structure of the organism (Polleux et al., 1997; Colombani et al., 2012; Tuazon and Mullins, 2015; McLaren and Steventon, 2021; Ebisuya and Briscoe, 2018). This coordination takes place across different chemo-mechanical scales and is locally encoded at the cellular level through genetic programmes. Also known as gene regulatory networks, these programmes confer the cell with the dynamical properties required for different cell differentiation processes (Delás and Briscoe, 2020; Reik, 2007). Notably, any inaccuracies in the timing of cellular dynamics can lead to critical errors that disrupt the healthy development of the embryo (Polleux et al., 1997; Colombani et al., 2012). Conversely, flexibility in the timing between tissues allows for the plasticity required for evolutionary traits to emerge. Given the crucial role of timing, a fundamental question arises in the field of developmental biology: how is the timing of cellular decisions precisely controlled?

One strategy to answer this question involves comparing homologous developmental processes across different species (McNamara, 2012), particularly during the phylotypic stage when embryos closely resemble each other (Duboule, 2003). These processes usually preserve the same function across species, manifested as the same sequential pattern of activation of gene expression. Recent research has shown that these developmental programmes, with evolutionarily conserved cis-regulatory architectures, manifest different species-specific timings (Halley, 2022; Lázaro et al., 2023; Linaro et al., 2019). These different timings can vary by an order of magnitude across distinct species, and can be observed in processes such as motor neuron differentiation (Rayon et al., 2020), the segmentation clock (Diaz-Cuadros et al., 2020; Matsuda et al., 2020), and corticogenesis (Iwata et al., 2023). During each one of these processes, the same genes are expressed in the same sequential manner, but at different speeds. Drawing an analogy from music, the species-specific execution speed of otherwise analogous genetic programmes has been termed the tempo of the system (Rayon et al., 2020; Diaz-Cuadros and Pourquié, 2023; Rayon, 2023; Simpson and Alberio, 2023).

Furthermore, the species-specific tempo in these processes can be inherently cell intrinsic, not requiring intercellular communication (Ebisuya and Briscoe, 2018; Busby and Steventon, 2021; Fuhrmann et al., 2020). This observation has resulted in an extensive search for the intracellular mechanisms governing tempo, with recent candidates focusing on specific molecular steps of gene regulation, such as protein stability, mRNA splicing, phosphorylation or ubiquitylation (Rayon et al., 2020; Matsuda et al., 2020; Blatnik et al., 2022; Sonnen, 2023; Iwata et al., 2023; Diaz-Cuadros et al., 2023). Contributing to this complexity, all these different mechanisms can affect each other, which, given the inherent non-linearity of gene regulation, impedes a straightforward mechanistic identification of tempo control.

The impact of non-linear feedback in gene regulation has been efficiently tackled in the past through mathematical modelling within the dynamical systems framework. The success of this strategy resides in its ability to connect molecular mechanisms with a general understanding of the compatible dynamics (Tomlin and Axelrod, 2007; Cooke and Zeeman, 1976; von Dassow et al., 2000; Dillon et al., 2003; Erban and Othmer, 2005). Yet the vast majority of the dynamical systems literature has been centred around steady-states, usually associated with available cell types, and very few studies have focused on the transient dynamics leading to those states (Tufcea and François, 2015; Zhu et al., 2017). This gap highlights the need for novel tools to comprehend the non-linear mechanisms governing tempo variations across species. Here, we devise a mathematical approach that connects the dynamic sequence of relative gene expression levels with the underlying biomolecular processes.

Before attempting to compare quantitatively the tempo between two species, it is necessary to ask whether they can actually be compared. Similarly, when we claim that the differentiation of two cells is the same except for their tempo, what do we mean by ‘the same’? In essence, we need to identify the differential characteristics of two temporal gene expression profiles that allow us to assert that they perform the same function but at a different tempo. Within the context of cell differentiation, one way to contain the function of a genetic programme is the conserved sequence of relative expression of different genes, independently of the speed at which this sequence is travelled. This definition is a useful proxy for biological function as it maintains many of the markers traditionally used to study embryonic development.

Adopting the jargon of dynamical systems to describe the situation in which two cells can have different gene expression trajectories, while still following the same orbits of gene expression, the orbits of the system can be visualised as the continuous sequence of cell states in the gene expression phase-space (Fig. 1A,B).

This identification of gene expression orbits with their core functional activity harmonises with the proposed concept of dynamical modularity proposed by Jaeger and Monk (2021): ‘A dynamical module or subsystem is defined by an elementary—distinguishable and quasi-autonomous—activity-function that corresponds to some dynamical regime within the broader dynamical repertoire of a complex system’. This way, the function of a system (or ‘dynamical module’) can be defined by ‘what the module does or, more precisely, by the particular kind of dynamics it generates’. It should be noted that other notions of function exist, such as the ‘use-function’, for which the function is determined by the downstream effects of the dynamical module, i.e. the effect that changes in the system have on the long-term development of the organism. We have decided against this definition in favour of activity-function, as we know from in vitro studies that tempo is cell intrinsic (unlike use-function), and because incorporating all the complexities of downstream effects would not be possible. Henceforth, here we will adopt Jaeger and Monk's ’activity-function’ definition when referring to ‘function’.

We can employ this definition to compare how close the functions of two different dynamical systems are by quantifying the similarity of their gene expression orbits (Fig. 1C). This strategy also reframes our exploration of tempo-controlling mechanisms: we can search for a suite of biochemical perturbations that sustain the integrity of a system's orbit. Thus, the effect of these biochemical perturbations can exclusively change the tempo. In addition, given the ability of mathematical models to incorporate different levels of biochemical detail in gene expression, we can use this framework to compare how different molecular mechanisms (e.g. post-translational modifications or promoter occupancy) affect the orbits of expression of given genes.

In general, it is not anticipated that two organisms will exhibit precisely identical orbits, as slight variations in orbits will maintain the overall function of the cell differentiation programme (Fig. 1C). Similarly, changes in biochemical rates are not expected to result in exact orbital equivalence (Eqn 2). Taken together, this highlights the requirement to define an orbital distance that measures the divergence between orbits in different systems.

Here, we will evaluate the distance of two orbits using the Fréchet distance, which captures the maximum deviation between two orbits, becoming zero when two systems are orbitally equivalent. The Fréchet distance between paths A and B can be intuitively understood as the following: the Fréchet distance is the minimum length of a leash required if a person were to walk along path A and their dog along path B, with both allowed to vary their speeds independently. A formal definition can be found in the supplementary Materials and Methods. Several alternative definitions of orbital distance are possible (Cleasby et al., 2019; Tao et al., 2021), each able to capture different aspects of the similarity between two different gene expression sequences (see supplementary Materials and Methods and Fig. S1 for more examples of possible measures). For example, a commonly used measure in biology is dynamical time warping (DTW). The main difference between Fréchet distance and DTW is that the latter is an aggregate method (accumulating differences along both orbits), whereas Fréchet is a bottleneck method, identifying the most extreme difference between two orbits. We focus on the former, because it prioritises penalties of local dramatic deformations of the orbit as opposed to a small orbit shift.

Ultimately, the choice of orbital distance needs to address how similar the functions of two systems are. This is not only required for two mathematical models to be compared, but also applies to the comparison of two experimental datasets of gene expression for which tempo is interrogated (regardless of whether the underlying gene regulatory network is known).

One immediate result from the orbital equivalence framework is that identical global rescaling of protein production and degradation rates of the transcription factors of a given gene regulatory network will preserve the orbit while controlling tempo. Interestingly, this rescaling does not have to be constant in time but can change arbitrarily in time (or even gene expression state). To illustrate this result, we analysed the dynamics of gene expression for a well-studied model of differentiation of neural progenitors in the ventral neural tube that includes the interaction of four genes (Fig. 3A) (Cohen et al., 2014; Rayon et al., 2020). Neural precursors initially express Pax6 and Irx3, and, when subjected to high Shh signal levels, undergo a transient gene expression cascade. This process involves the orderly activation/deactivation of genes within the network, culminating in the expression of Nkx2-2, which steers the precursors towards a V3 interneuron identity. We modified the model by controlling a prefactor, μ(t), multiplying the production and degradation rates (Fig. 3B). As expected, different temporal profiles of μ(t) yield different trajectories. Nevertheless, the orbits of gene expression are conserved, i.e. the sequence of relative expression of all the genes involved in the neuron differentiation process is unaltered, including the final expression state.

Hence, this rescaling provides a general mechanism to preserve the orbit of gene expression independently of changes in rates along the cell cycle or through external perturbations as long as production and degradation are rescaled in the same way. This result applies to any given regulatory network independently of the complexity of its topology, suggesting a robust mechanism to tune tempo of a regulatory system, allowing for evolutionary strategies.

Recent experimental work by Rayon et al. (2020) on neural precursor differentiation used the same model to demonstrate that protein degradation rate can control tempo. In Fig. 3C, we compare the time factor resulting from a modification of degradation rates only (as was done in the Rayon et al. study), versus a global rescaling of production and degradation. As expected, we observed that global rescaling affects the tempo linearly, and indefinitely, achieving any possible tempo. By contrast, modifying degradation rates alone restricts the range of achievable tempos, meaning that the Nkx2-2 state cannot be achieved when the change in degradation rate exceeds a twofold increase. This failure stems from the fact that isolated changes in protein degradation rate alone cannot maintain the shape of the orbit (Fig. 3D). In this case, it can even change the final steady-state of the trajectory. This phenomenon is further observed by an increase in orbital distance, which grows as the degradation rate changes. A significant jump in the orbital distance occurs when the Nkx2-2 state becomes unattainable, thereby altering the cell fate to a high Olig2 (motor neuron) fate.

Therefore, although changing degradation rates can explain the time factor difference of 2.5 between human and mouse reported by Rayon et al. (2020), it cannot immediately explain the greater differences present between other species. For example, zebrafish neural precursor differentiation lasts less than a day, whereas for humans it takes around two weeks (Rayon and Briscoe, 2021), a time factor of 10 or more. Supporting our predicted general mechanism of tempo control, recent experimental studies investigating the species-specific timing of the segmentation clock (Matsuda et al., 2020) have identified differences in protein stability and production rates across species. Although we cannot conclude that there are other molecular mechanisms at play allowing the extension of tempo beyond the limits allowed by changes in degradation rates, our framework proposes a straightforward molecular strategy.

In general, arbitrary perturbations other than global scalings of production and degradation will not necessarily result in orbital equivalence. To showcase the applicability of the orbital equivalence framework in these scenarios, we implemented it within a well-known gene regulatory network, the repressilator (Fig. 4A). Despite not corresponding with a particular developmental biology example, it is a paradigmatic circuit capable of gene expression oscillations that offer the flexibility to incorporate different degrees of biochemical complexity (Perez-Carrasco et al., 2018; Elowitz and Leibler, 2000; Wu et al., 2017) and different tempos (Zhang et al., 2022). In addition, sustained oscillations result in closed orbits in the state-space (Fig. 4B,C), allowing for a clearer visualisation of the concept of orbital distance.

To explore the dynamical properties available through changes in the directionality of reactions towards the dimer state r_i, we undertook a parameter search for different perturbations of the parameters r_i with respect to an arbitrarily chosen reference system. Analysis of the resulting oscillatory systems shows high variability in the possible oscillatory behaviours (Fig. 4D,E). To visualise the results of this search, we can pick out a subset of systems from Fig. 4E (highlighted in orange) and display their orbits (Fig. 4F) and trajectories (Fig. 4G) in a grid reflecting their relative position in Fig. 4E. For both Fig. 4F and G, rows correspond to systems with the same period (tempo), and columns correspond to systems with the same orbital distance to the reference system, which increases as you go from left to right. The resulting orbits not only exhibit a large variability of shapes (Fig. 4F), but there is also variability in the periods of the orbits (Fig. 4G). Most interestingly, there are many cases in which the period of the orbit can be perturbed without changing the orbit of the system (see the left-most column of Fig. 4F,G). Similarly, there are cases in which different orbits are obtained while preserving the period of the oscillations, suggesting that there are orthogonal molecular strategies to change orbit shape and tempo.

Observing the variability of behaviours accessible by changes in r_i, we investigated whether we could predict the orbital distance and resulting period of the perturbed systems obtained in our parameter search by analytical inspection of the prefactors μ_i (Eqn 11). Using our analytical definition of prefactor heterogeneity D_σ (as outlined in Eqn 4 and Eqn 5) and comparing it with the orbital distance calculated by numerical integrations of the perturbed systems, we observed a linear relationship between the two (Fig. 5A). This relationship reveals an analytical means to predict the magnitude of the changes in the orbit's shape induced by a perturbation in the repressilator, without needing to integrate the perturbed ordinary differential equations.

Similarly, we investigated whether the period of oscillations in the perturbed system could be predicted by analysing only the prefactor of the new system. In this analysis, we found that the average prefactor serves as a reliable predictor of oscillation periods (Fig. 5B). This finding offers a method for forecasting the resulting oscillatory period without the numerical necessity of integrating the ordinary differential equations of the perturbed system.

Although the μ_i prefactors are good indicators of the effect of a perturbation in the system, they come with two main disadvantages. The prefactor needs to be evaluated along the orbit of the reference system (Eqn 5), and the prefactor is not immediately linked with the molecular mechanism in charge of orbital distance and tempo. Hence, we investigated whether we can identify molecular mechanisms of tempo and orbit shape control with definitions that link directly to the combined rates towards dimerisation r_i. Initial inspection of the role that individual rates have on the prefactor reveals that, although there is some correlation between each r_i and the resulting prefactor heterogeneity and average prefactor, there is not a strong predictive relationship (Fig. 5C,D). These results suggest that the effect that changes of prefactor have on the orbit originate from a combination of perturbations in the ensemble of rates included in each r_i. In particular, we identified that the average of the composed dimerisation rates can be used as a good predictor of the average prefactors, and consequently the period of the resulting orbit (Fig. 5E, Fig. S2B). However, we also find that the heterogeneity of r_i allows us to predict the prefactor heterogeneity, resulting in a predictor of the orbital distance (Fig. 5E, Fig. S2C). Therefore, our framework has allowed us to find underlying mechanisms (r₁, r₂ and r₃) that affect orbital distance. Unlike the orbital distance, or prefactor analysis (e.g. D_σ), calculating r_i heterogeneity and average does not require any knowledge of the orbit of the reference system. Instead, it provides direct molecular mechanistic insight into the effects of post-translational rates on trajectories.

These results provide two orthogonal molecular mechanisms to control orbit shape and tempo. Hence, this allows us to use these mechanisms to design sets of parameters that will preserve the orbit of a system while exclusively changing the tempo of the system, independently of the reference system. We can achieve this by tuning the average r_i while fixing the heterogeneity in the prefactor. Implementing this strategy, we successfully gained the capability to fine-tune the tempo for orbits of various shapes (Fig. 5F), resulting in achievable speed differences spanning a full order of magnitude.

We have developed a mathematical framework using dynamical systems that can be used for investigating molecular mechanisms in charge of inter-species differences in developmental tempo. This framework is based on the notion of orbital distance, assigning a dynamical function to a system based on a defined progression of relative gene expression states. We have shown how we can explore this idea in two parallel ways: (1) by comparing orbits directly from gene expression data, and (2) by providing an analytical toolkit connecting ordinary differential equation representations with developmental tempo.

The versatility of the techniques introduced means that they can be applied to many different genetic programmes in development. The similarity measures based on gene expression data can be immediately used to compare function between organisms regardless of the underlying gene regulatory network. Although this paper employs the Fréchet distance for this comparison, many other possible distances are possible (see supplementary Materials and Methods), reflecting different definitions of function in embryogenesis. Similarly, we focus our definition of ‘function’ as the activity-function resulting from the orbit of a dynamical system, but identical orbits with different tempo can still entail different phenotypes (e.g. patterns of different width in the oscillatory example of Fig. 1), requiring the concept of function to be revisited when employing the framework in different systems.

The mathematical framework becomes very useful when contrasted with more traditional representations of the differentiation process as trajectories in a Waddingtonian landscape. Recent advances in this field have focused on reconstructing these landscapes from experimental data, usually fixing the location and change of steady states and basins of attraction of the landscape (Sáez et al., 2022; Verd et al., 2014). Here, we show that there are infinite available transformations of the dynamical system that preserve the same identical landscape. This transformation can be achieved through a prefactor of the dynamics (Fig. 2A). Overall, this reveals two compatible formalisms that separate the geometry and the tempo of a dynamical system's landscape. This not only allows us to understand the underlying mechanisms of tempo control, but also informs how we interpret and find landscape representations of a particular differentiation process.

Although our concept of orbital equivalence focuses on a preservation of the direction of the flow along a typical differentiation trajectory, further research should be done to include the robustness of these orbits under noise, delays, or changes in the initial conditions. In this latter case, the concept of orbital equivalence is closely related to the concept of homeorhesis (Verd et al., 2014; Young et al., 2017; Tufcea and François, 2015), whereby the relevant orbit of the system is described by the slow manifold of the dynamical system, buffering the effect that changes in the initial conditions have on the resulting orbit.

One of the main features of this study is showcasing that tempo can be controlled through intermediate reaction steps during the regulatory process. Previous experimental studies have shown that inter-species developmental tempo depends heavily on specific biochemical parameters (Matsuda et al., 2020; Rayon et al., 2020; Lázaro et al., 2023; Filina et al., 2022). Thus, the tools developed could be employed to investigate such systems to understand the impact and limits that different parameter perturbations will have on the resulting trajectories of gene expression. In order to do so, quasi-steady-state approximations need to be developed that can capture the impact that metabolic and other biochemical processes have on the overall change of gene expression in the state-space. In our case, we identified dimerisation rates as a pathway to control tempo; nevertheless, we anticipate that other intermediate steps may have similar temporal implications.

All in all, this framework not only shows a route to quantitatively connect molecular mechanisms with tempo in development, but it can also be extended to other areas. For instance, in an evolutionary context, the same tools can be used to identify possible molecular routes to change tempo while preserving the function of a system, hence identifying evolutionary strategies that showcase the molecular requirements at different tempos. Similarly, these strategies can be employed directly in synthetic circuit design, with the aim of building circuits that allow for tempo control of specific functions, or even controlling the tempo of in vitro cell differentiation.

We thank Dr Smitha Maretvadakethope for her feedback on the manuscript. Our acknowledgements also extend to Prof. Nick Jones and Prof. Robert Endres for their insightful discussions about the project. R.P.-C. extends gratitude for the discussions at the KITP, which were supported in part by the National Science Foundation (NSF PHY-1748958 and PHY-2309135).

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.202950.reviewer-comments.pdf