In our paper, ‘A theoretical framework for the regulation of Shh morphogen-controlled gene expression’ by Cohen et al. (2014) we formulate a mathematical model of gene regulation by morphogen signalling that brings together empirical findings from several sources including Balaskas et al. (2012), Oosterveen et al. (2012, 2013) and Peterson et al. (2012). We use an approach based on statistical thermodynamic ensemble models of gene regulation and Approximate Bayesian Computation. We argue that the mathematical model provides a single coherent framework that explains experimental observations and that the approach can be applied to similar morphogen systems.
Uhde and Ericson (2016) do not dispute our mathematical model. Instead they claim: (1) our study ‘overlooks previous experimental work’ and we have not ‘acknowledged or discussed…conceptual conclusions’ of the Ericson lab; and (2) our conclusions are ‘conceptually indistinguishable’ from Oosterveen et al. (2013, 2012) and Peterson et al. (2012), and our ‘major conclusions primarily confirm published models of Shh interpretation’.
We disagree. Here, we clarify the issues that appear to have caused these misunderstandings.
The contribution of the Ericson and McMahon laboratories to understanding the roles and dissecting the architecture of key patterning genes in the neural tube is indisputable. Far from overlooking their studies, we use these explicitly as motivation for the mathematical model and we cite the papers extensively throughout the manuscript (we cite the Peterson study seven times and the Oosterveen papers nine times). We do not think readers could fail to notice our references to these papers or easily miss our use of their work to support the model. Moreover, as Uhde and Ericson (2016) point out, we also reviewed and cited their work extensively in a separate publication (Cohen et al., 2013).
Nevertheless, the description of morphogen-controlled gene expression is the product of a large number of studies that go well beyond the neural tube field (Briscoe and Small, 2015). The model of gene regulation we construct contains three types of input: (1) uniformly expressed transcription factors (TFs). The role of these in morphogen interpretation became obvious from work on Zelda and Stat92E in Drosophila (e.g.Kanodia et al., 2012; Nien et al., 2011). (2) A transcriptional effector of the morphogen for which there is a lack of correlation between binding affinity and the position of target gene activation. This has been extensively documented for Bicoid in the Gap gene system (Ochoa-Espinosa et al., 2005). (3) A set of morphogen-regulated TFs that form a transcriptional network. The Gap genes also provide a well-established example of the importance of transcriptional network dynamics in morphogen pattern formation (e.g. Jaeger et al., 2004; Manu et al., 2009).
Each of these elements was described in Drosophila prior to the work of Oosterveen et al. (2012, 2013) and Peterson et al. (2012). Furthermore, the idea that cis-regulatory modules combine multiple inputs to ‘compute’ an output is at the heart of the gene regulatory network framework developed by Davidson and colleagues (reviewed in Davidson, 2010). Thus the precedents for the broad conceptual conclusions to which Uhde and Ericson (2016) refer arose from studies of non-vertebrate systems that predate the work of Oosterveen et al. (2012, 2013) and Peterson et al. (2012).
Our intention by citing the non-vertebrate studies was not to diminish the contribution of the Ericson and McMahon labs, but to provide a broader context and the appropriate background. For example, in the section dealing with the function of uniformly expressed TFs we write: ‘…previous studies have demonstrated how the levels of binding of a spatially uniform factor to target genes in a morphogen patterning system can significantly influence their expression profiles (Kanodia et al., 2012). In the neural tube, the TF Sox2 has been suggested to provide a spatially uniform activation input into neurally expressed genes (Bailey et al., 2006; Oosterveen et al., 2012; Peterson et al., 2012)’. We believe that these comparisons are highly relevant and emphasize the importance of Oosterveen et al. (2012, 2013) and Peterson et al. (2012). Taken together, the studies suggest common principles underpin the transcriptional interpretation of morphogen signalling in several tissues.
More importantly, the suggestion that our conclusions are ‘conceptually indistinguishable’ and that our ‘major conclusions primarily confirm published models of Shh interpretation’ misses the key point of our paper. Cohen et al. (2014) describe and analyse a mathematical model. The Oosterveen and Peterson studies do not contain mathematical models, neither does Cohen et al. (2013). Moreover, the interpretation that Oosterveen et al. (2012) offer of their data is not equivalent to the mathematical model in Cohen et al. (2014).
We believe that the explanatory and predictive power of mathematical models is of most value when firmly rooted in experimental observations. The empirical observations we use to construct the model are based on the studies of Oosterveen et al. (2013, 2012), Peterson et al. (2012), among many other studies, and we cite these papers accordingly. These observations are the basis for the model not, as Uhde and Ericson (2016) seem to suggest, the ‘conclusions’ of the model. In our view, the analytical framework in Cohen et al. (2014) helps to rigorously establish the relationships between different pieces of empirical evidence and formulates a mechanistic, predictive model of gene regulation.
Moreover, mathematical models, particularly of non-linear dynamical systems such as transcriptional networks, often provide insight into complex behaviours that are difficult to discern from experiment alone. We think this is the case here. As a consequence, there are several differences that distinguish the model proposed in Cohen et al. (2014) from Oosterveen et al. (2012, 2013).
(1) In the Cohen et al. (2014) model there are no inherent differences between how target genes interpret GliA and GliR gradient. By contrast, the study by Oosterveen et al. (2012) proposes two classes of genes: ‘local’ genes that interpret ‘the balance between GliA and GliR’, and ‘long-range’ genes that only interpret GliR and have more dorsally positioned boundaries. They label this a ‘GliR gradient interpretation model’ and comment: ‘This GliR-gradient interpretation model differs significantly from prevailing models suggesting…that cells strictly measure the balance between GliA and GliR’. This leads Oosterveen et al. (2012) to their principal conclusion that the interpretation of Shh signalling involves ‘mechanistic differences’ between the ‘local’ and ‘long range’.
In the model proposed in Cohen et al. (2014), the expression of both short and long-range Gli-regulated genes depend on the concentration and strength of binding of both GliR and GliA (the binding affinity of both is the same – parameterized by a single value). This could be summarized as ‘genes measure the balance between GliA and GliR’. Thus, Cohen et al. (2014) do not invoke two classes of genes nor mechanistic differences between short and long-range target genes.
The model proposed in Cohen et al. (2014) also reveals that differences in GBS affinity between target genes are not necessary to explain the observed spatial-temporal dynamics of target gene expression. This appears at odds with the proposal of ‘mechanistic differences’ in the transcriptional regulation of ‘local’ and ‘long-range’ target genes (Oosterveen et al., 2012).
(2) Uhde and Ericson (2016) state that their data indicate that ‘Gli-mediated gene activation [is] a largely concentration-independent event’. This is not the case in the mathematical model described in Cohen et al. (2014) in which the response of target genes is dependent on the concentrations of GliA and GliR. The way a specific target gene responds to alterations in activator and repressor levels depends not only on the Gli input but also its other inputs. An important consequence of this is that changes in the Gli binding affinity of a target gene can have opposite effects on the range of activation of different genes. We believe that this represents a good example of how the mathematical model helps explain the experimental data. We make this point in our manuscript by citing data in Oosterveen et al. (2012) and Peterson et al. (2012).
(3) As Uhde and Ericson (2016) indicate, our model describes a ‘neutral point’. However, our definition of the neutral point appears to differ conceptually from their interpretation. In their correspondence, Uhde and Ericson imply that the neutral point is a specific location in the tissue between the Nkx2.2 and Olig2 boundaries. Furthermore, they suggest that genes on either side of this point have different Gli binding affinities. In the mathematical model this is not the case. In Cohen et al. (2014), we define the neutral point for a gene as the point in the GliA/GliR gradient at which altering the GBS affinity does not alter the probability of gene expression (see eqn 5). We show how this point depends on the concentration of GliA/GliR and the basal levels of expression of each gene (see eqn 6 and eqn S4). Hence the neutral point is independent of the binding affinity for Gli TFs (parameter K) and is not a single position in the tissue – for each target gene it depends on the non-Gli regulatory input (eqn S4). Importantly, our analysis indicates that the observation that changes in GBS affinity result in opposite shifts of gene expression boundaries on either side of the neutral point is an emergent property of the model.
We agree with Uhde and Ericson that this behaviour is evident in the Oosterveen and Peterson experimental data, but we could not find the idea of a ‘neutral point’ proposed in Oosterveen et al. (2012, 2013). Instead Oosterveen et al. (2012) use the data to propose ‘mechanistic differences’ between ‘local’ and ‘long-range’ interpretations of Shh signalling and that ‘Gli activators have a noninstructive role’. In Cohen et al. (2014), we highlight and cite the experimental evidence of the neutral point revealed by their experimental data and show how a mathematical model suggests a single mechanism to explain the experimental observations. We think this represents one of the successes of the mathematical model and it illustrates how such models provide new insight into experimental data.
(4) Oosterveen et al. (2012) suggest that gene regulation involves ‘cooperative’ interactions between Gli and HD proteins and between Gli and SoxB1. The model we formulate does not contain these cooperative interactions. This does not rule out cooperative interactions in vivo. However, these are not required in the mathematical model.
(5) Uhde and Ericson (2016) say that we ‘incorrectly state “…analysis of GBSs within enhancers of Shh target genes failed to find a positive correlation between binding site affinity and range of gene induction (Oosterveen et al., 2012; Peterson et al., 2012)” ’. We are confused by this statement. The data reported in fig. 3A of Oosterveen et al. (2012) and table 1 of Peterson et al. (2012) lack a positive correlation between binding site affinity and range of gene induction – at best there might be a weak negative correlation in the Oosterveen data, which is not evident in the Peterson analysis (see Fig. 1). In Oosterveen et al. (2012) two classes of genes, ‘local’ and ‘long range’ are defined, but even within each of these two classes Oosterveen et al. (2012) conclude ‘there is no predictive correlation between gene expression pattern and affinity score or number of GBSs’ (p. 1009). This appears to be in line with our statement. We note that the P19 data in Oosterveen et al. (2012) also support this conclusion – fig. 3D in Oosterveen et al. (2012) shows no difference between the Nkx2.2 and FoxA2 activity and the statistical significance of the other differences is unclear.
Scatter plots of the predicted binding affinities of the putative gli binding sites (GBSs) associated with the indicated genes as reported in Oosterveen et al. (2012) and Peterson et al. (2012) . Data were extracted from fig. 3A of Oosterveen et al. (2012) and from table 1 of Peterson et al. (2012).
Scatter plots of the predicted binding affinities of the putative gli binding sites (GBSs) associated with the indicated genes as reported in Oosterveen et al. (2012) and Peterson et al. (2012) . Data were extracted from fig. 3A of Oosterveen et al. (2012) and from table 1 of Peterson et al. (2012).
While not a point of distinction between the conclusions of Cohen et al. and the earlier studies, Uhde and Ericson (2016) raise concerns about our discussion of hysteresis. The term ‘hysteresis' was originally coined to describe the behaviour of magnetic materials and is widely used in dynamical systems theory to indicate that the output of a system depends not only on its current input, but also on past inputs. We show that the mathematical model developed in Cohen et al. (2014) displays similar hysteresis to the simpler dynamical system described in Balaskas et al. (2012). This is a validation of the model formulism and it links the Cohen et al. model to the experiments performed in Balaskas et al. (2012). Hysteresis suggests an explanation for the maintenance of gene expression in cells in which the levels of Gli activity decrease over time; an observation described in Balaskas et al. (2012) and recently independently observed in the elegant quantitative approach taken by Junker et al. (2014). This mechanism is distinct from that proposed in Lek et al. (2010), which is why we did not cite that work.
In summary, there are several features that distinguish the mathematical model described in Cohen et al. (2014) from previous interpretations. The Cohen et al. (2014) model demonstrates that neither differences in the binding affinity of GliA and GliR, nor differences in GBS affinity between target genes are required to explain the patterns of gene expression. Contrary to the conclusions of Oosterveen et al. (2012), the model in Cohen et al. (2014) proposes that all target genes respond to the ratio of GliA and GliR and the mathematical model does not impose mechanistic differences between local and long-range targets. Together, the Cohen et al. model provides a distinct mechanistic explanation for the experimentally observed position-dependent shifts in gene expression upon perturbations of binding site affinity.
We note that the section entitled ‘Discussion’ in Cohen et al. (2014) was originally titled ‘Conclusion’ to signify it as a short summary and the section entitled ‘Results’ was originally titled ‘Results and Discussion’. These headings were changed in response to an editorial request, after acceptance, to match Development's house style and we regret any confusion this caused.
We are saddened that Uhde and Ericson believe we undervalue their work. The intention of Cohen et al. (2014) was not to diminish or undermine studies to which they and others have contributed. Rather, we believe that accommodating the empirical data of Oosterveen et al. (2012, 2013) and Peterson et al. (2012) in a single theoretical framework and reconciling this with other studies in the field emphasizes the importance and success of their work. We think this exchange of correspondence also highlights the benefits of mathematical models: they provide formal, precise and transparent descriptions of ideas that are not subject to the ambiguities or differences in interpretations of narrative accounts or informal ‘cartoon’ models. Moreover, mathematical models make clear predictions that can be tested experimentally and we look forward to working with the Ericson lab and others in the field to revise, extend or refute current models.
