Autoregulation of Transcription and Translation: A Qualitative Analysis

The regulation of both mRNA transcription and translation by down-stream gene products allows for a range of rich dynamical behaviours (e.g. homeostatic, oscillatory, excitability and intermittent solutions). Here, qualitative analysis is applied to an existing model of a gene regulatory network in which a protein dimer inhibits its own transcription and upregulates its own translation rate. It is demonstrated that the model possesses a unique steady state, conditions are derived under which limit cycle solutions arise and estimates are provided for the oscillator period in the limiting case of a relaxation oscillator. The analysis demonstrates that oscillations can arise only if mRNA is more stable than protein and the effect of nonlinear translation inhibition is sufficiently strong. Moreover, it is shown that the oscillation period can vary non-monotonically with transcription rate. Thus the proposed framework can provide an explanation for observed species-specific dependency of segmentation clock period on Notch signalling activity. Finally, this study facilitates the application of the proposed model to more general biological settings where post transcriptional regulation effects are likely important.

A conserved principle of the Hes/Her oscillator, now known to be present in many different cell types (Kageyama et al. 2007), is that dimerised members of the basic Helix-loop-helix family of transcription factors (e.g. Hes7, Hes1, Her7) inhibit their own transcription and therefore provide a negative feedback loop. The Notch signalling pathway, which plays a crucial role in embryo development, tissue homeostasis (van Es et al. 2005) and cancer (Mollen et al. 2018;Allenspach et al. 2002;Siebel and Lendahl 2017), can activate the transcription of Hes/Her genes. During canonical in trans Notch signalling, a Notch ligand in a signalling cell activates a Notch receptor in a neighbour, resulting in the release of the Notch intracellular domain (NICD) in the receiver, which regulates the transcription of Notch target genes. As at least in some biological contexts, such as the segmentation clock, Notch receptors are themselves a target of Notch signalling and levels of the Delta ligand can be regulated by Hes7 (Bone et al. 2014), the study of Notch signalling is a highly nonlinear problem.
Upon inclusion of time delays that represent processes such as transcription, splicing, transport and translation, it has been shown that negative feedback of transcription is sufficient to give rise to oscillations (Lewis 2003;Monk 2003). Moreover, it has been shown that the spatial diffusion of mRNA and protein is a sufficient mechanism to give rise to oscillations in a negative feedback system (Sturrock et al. 2011;Chaplain et al. 2015). Each of the above models makes the assumption that the translation of mRNA is linear, and thus unregulated.
Recent experimental observations challenge, at least in specific biological contexts, many existing models of the Notch signalling pathway. Oates and coworkers have demonstrated that when levels of Delta ligand are increased in presomitic mesoderm (PSM) cells, the tissue scale oscillator period decreases (Liao et al. 2016). Moreover, when levels of Notch signalling are reduced via treatment with the gamma secretase inhibitor DAPT, which blocks the release of NICD, the tissue scale oscillator period increases (Herrgen et al. 2010). Thus in the zebrafish embryo, the tissue-scale oscillator period appears to be anticorrelated with Notch signalling activity. In contrast, Dale and coworkers have demonstrated that when mouse and chick embryos are exposed to pharmacological treatments that increase levels of NICD, the tissue scale period increases (Wiedermann et al. 2015). Notably, a prediction of the delayed feedback models of the Her oscillator (Lewis 2003) is that the clock period has a strong dependence on the mRNA and protein half lives and time delays but not on transcription or translation rates (Lewis 2003).
Suggestions that mouse PSM tissue behaves like an excitable medium are also difficult to reconcile with delayed negative feedback models of the Notch signalling pathway. It has been identified that NICD is necessary for the oscillations of the segmentation clock in the presence of mechanosensitive Yap signalling (Hubaud et al. 2017). However, when Yap signalling is pharmacologically inhibited, oscillations could still proceed in the absence of Notch signalling. The presence of a Yap-signalling dependent threshold led the authors to conclude that the system under study behaved like an excitable medium. However, there is currently no molecular scale model of Hes7 dynamics that can account for such excitability.
In the Hes1 oscillator in mouse neural cells, the miRNA mir-9 has been identified (Bonev et al. 2012;Goodfellow et al. 2014) as a part of a double negative (i.e. positive) feedback loop in which mir-9 is under the same transcriptional control as the Hes1 gene but serves to inhibit translation. Together, these observations indicate that, at least in some specific biological contexts, the negative feedback model of the Notch signalling pathway is incomplete.
We recently developed an ordinary differential equation model that postulates that an intermediary, X , that is under the same transcriptional control as mRNA, M, inhibits translation [see Fig. 1a (Murray et al. 2021)]. Assuming quasi equilibrium for X , these assumptions introduce a positive feedback loop such that the translation rate increases sigmoidally with protein concentration (see Fig. 1b). Letting M = M(t) and P = P(t) represent the concentrations of mRNA and its corresponding protein at time t, respectively, the governing ODEs are where k 1 is the maximal transcription rate, P 0 is the protein concentration at which transcription rate is half maximal, k 2 is the mRNA degradation rate, k 3 is the basal translation rate, k 4 is a translation rate that is inhibited by X , α is maximal level of X at steady state, X 0 is an IC50 constant for translational activation and k 5 is the protein degradation rate.
Using parameter values based on the zebrafish Her oscillator, it was shown, using numerical exploration, that Eq. (1) can possess excitable, homeostatic or oscillatory solutions (Murray et al. 2021). Using numerical continuation it was shown that Eq. (1) possess a subcritical Hopf bifurcation such that, in a particular region of parameter space, unstable limit cycle, stable limit cycle and stable steady state solutions coexist. In this case a stochastic implementation of the model is capable of exhibiting intermittent oscillations whereby noise switches the dynamics between a stable limit cycle and a stable steady state. Finally, it was shown that the oscillator period has, for the considered parameters, an inverse dependence on the transcription rate k 1 . Hence the proposal that regulation of translation, as well as transcription, rates provides a minimal framework that yields phenomena consistent with recent experimental observations. Whilst the previous work used numerical solutions to demonstrate interesting model behaviours, a qualitative analysis of the model behaviour is required in order that the model can be explored in more general biological contexts. Here this issue is addressed. The approach taken allows one to relate different Notch signalling behaviours in species-specific contexts (e.g. in which reaction rates may differ significantly). Parameter regimes are identified in which one expects to find different modes of behaviour Translation and transcription rates plotted against protein levels, P (e.g. excitability, homeostasis, oscillations). Finally, an estimate is derived for the oscillator period and amplitude in the relaxation oscillator limit.

Nondimensionalisation
Consider the dimensionless variables where Note that time has been nondimensionalised on the protein degradation timescale, the parameter η 4 represents the strength of the sigmoidal effect on translation rate and that η 3 is the ratio of mRNA to protein degradation rates. See Table 1 for typical values.

Nullclines
The p nullcline, given bym has two distinct real positive turning points if the condition η 4 < 1 9 holds (see Appendix A). The turning points occur at approximately and ( p 2 , m 2 ) = (1, 2), (see Fig. 2). Note that the condition η 4 < 1/9 implies that m 1 > m 2 and p 1 < p 2 . Hence ( p 1 , m 1 ) is a local maximum and ( p 2 , m 2 ) is a local minimum. The m nullcline, given bym is monotonically decreasing for p > 0 with an IC50 at p = √ η 2 and local maximum of η 1 /η 3 . In order that nondimensional parameters correspond to positive dimensional parameters, the condition η 2 < η 4  Table 1 must hold (see Appendix A). Hence the IC50 for transcriptional inhibition must be at least an order of magnitude less than the IC50 for the translational switch (which occurs at approximately p = 1).

Steady State Analysis
Suppose that (m * , p * ) is a steady state of Eq. (2). Upon elimination of m * , p * satisfies the fifth order polynomial Recalling that η j > 0 ∀ j, application of Descartes' rule of signs implies that there are at most three real positive solutions of Eq. (6). Moreover, as Equation (6)

Linear Stability Analysis
After substitution for the identity the Jacobian matrix of equations (2) takes the form Given that Eq.
(3) possesses a unique steady state in the positive quadrant, the sign structure of the Jacobian matrix is given by

Intersections on the Left and Right Branches are Linearly Stable
Negativity of the (2, 2) entry of the Jacobian matrix implies that the left-hand side of which has previously been used to compute the turning points of the p nullcline (labelled as p 1 and p 2 , see Eq. (19) in Appendix A). Hence when p * < p 1 or p * > p 2 , such that the intersection between the m and p nullclines occurs on the left-or right-most branches of the p nullcline, respectively, the (2, 2) entry of the Jacobian matrix is negative and the steady state of Eq.

The Steady State on the Central Branch of the p Nullcline is Conditionally Linearly Stable
The determinant of the Jacobian matrix is positive definite (see Appendix C). Hence the unique steady state of Eq.
(2) is linearly unstable if and only if This inequality can be expressed as with the boundaries of the solution interval given by For a real and positive solution interval it is therefore required that which can, upon rearrangement, be written as Considering the case where η 3 < 1, a necessary (but not sufficient) condition for instability of the steady state is In summary, when the conditions hold, there is always a real interval of p * within which tr(A) is positive and the steady state is therefore linearly unstable. As p * is a monotonically increasing function of η 1 (see Appendix B) a corresponding interval of the parameter η 1 can always be found such that the unique steady state ( p * , m * ) is linearly unstable. This result implies that too little or too much basal transcription (i.e. k 1 ) will result in the disappearance of oscillatory solutions.

Limit Cycle Solutions
A confined set can be defined for Eq.
(2) (see Appendix D). Given the existence of a unique steady state, the Poincare Bendixson theorem can be applied in order to show that there is an interval of the parameter η 1 for which Eq.
(2) have limit cycle solutions.  Table 1 unless otherwise stated (Color figure online)

Numerical Continuation
Numerical continuation was performed (Dankowicz and Schilder 2013) to confirm the presence of a family of Hopf bifurcations in the η 1 − η 3 plane (see Fig. 3a). These numerical results indicate that, given that inequalities (14) hold, one can find an interval of the parameter η 1 in which there are limit cycle solutions. Note that the derived upper bound on η 3 , given by inequality (13), is consistent with the upper bound estimated using continuation. Moreover, by imposing the conditions such that the nullclines intersect in the middle branch of the p nullcline, a necessary condition for limit cycle solutions is These bounds are presented in Fig. 3a. In Fig. 3b, c the Hopf bifurcation surface is projected onto the η 1 − η 3 plane for different values of η 4 and η 2 , respectively. Note that, as expected, the maximum value of η 3 for which oscillatory solutions are possible varies with the parameter η 4 (see Fig. 3b) but not η 2 (see Fig. 3c). Numerical continuation also indicates that the classification of the Hopf bifurcation that arises for smaller η 1 is dependent on the parameter η 3 . For larger η 3 , there are two supercritical Hopf bifurcations. Here the amplitude of oscillations increases close to both bifurcation points (see Fig. 4a-c). However, for smaller η 3 the Hopf bifurcation is subcritical and one observes the emergence of a saddle node bifurcation of the limit cycle. In this case there is an interval of the parameter η 1 in which there is an unstable limit cycle, a stable steady state and a stable limit cycle (see Fig. 4d-f). In Bottom row, η 3 = 0.02. Left column, steady state levels of protein, p, plotted against η 1 (blue dashed line, unstable; red line, stable). Solid black lines represent maxima/minima of the p component of limit cycle solutions. Middle column-inset for left column. Right column-oscillator period is plotted against η 1 . Red markers-stable limit cycle. Blue markers-unstable limit cycle. Parameters as in Table 1 unless otherwise stated (Color figure online) the limiting case, where both η 1 and η 3 are small, the time scale of mRNA production and degradation are relatively long and the system behaves like a relaxation oscillator (see Fig. 4g-i). Notably, the dependence of the oscillator period on the parameter η 1 is in general not monotonic.

Period and Amplitude Estimate in the Relaxation Oscillator Limit
Under the assumption that inequalities (14) hold, an estimate for the oscillator period can be derived. Consider the case where the time scale of mRNA transcription and degradation is much longer than that of translation. After applying a fast-slow time scale analysis, where the mRNA is the slow variable, the limit cycle is approximated by a trajectory ABCD (see Fig. 5b) with coordinates A lower bound for the oscillator period (see Appendix E) is given by Thus in the relaxation oscillator limit the period varies inversely with the parameter η 1 . The amplitudes of protein and mRNA oscillation are approximated by respectively. In Fig. 5c-f the derived estimates for the oscillator period are compared with numerical estimates. It is noted that as the oscillator is made less stiff, a correction is needed to Eq. (17) that accounts for time spent close to the local maximum of the p nullcline (see Appendix E). In this case the estimate of the oscillator period no longer depends monotonically on the parameter η 1 .
Thus for the p nullcline to have two turning points there must be a significant upregulation of the net translation rate and the maximal level of X must be much larger than the IC50 for the upregulation of the translation rate. Upon expansion in the small  Table 1 (Color figure online) parameter η 4 , condition (13) can be approximated by These conditions imply the more restrictive bounds α X 0 > 16 and k 4 k 3 > 16.
A region of parameter space in which oscillations are possible is depicted in Fig. 6. The results imply that an experimental perturbation that independently either: (i) decreases the translation rate ratio; (ii) decreases the steady state level of X; or (iii) decreases mRNA stability relative to protein stability could be sufficient to move the system out of the oscillatory regime. Moreover, mRNA must be more stable the protein in order for oscillatory solutions to be possible. Upon redimensionalising Eq. (17), an estimate for the oscillator period in the relaxation oscillator limit is given by Notably, whilst the oscillator period increases linearly with the mRNA half-life (ln 2/k 2 ), it has an inverse dependence on the protein half life (ln 2/k 5 ). It is also inversely dependent on the transcription rate, k 1 , and there is a strong nonlinear dependence on translation rates (k 3 and k 4 ). The dimensional protein and mRNA oscillator amplitudes are given at leading order in the relaxation oscillator limit by respectively. Notably, the amplitude of protein and mRNA oscillations are independent of mRNA production and degradation rates. Moreover, the ratio of protein to mRNA amplitudes can be approximated by

Discussion
A model of the Notch signalling pathway was recently developed in which it was assumed that an intermediate factor that is under the same transcriptional regulation as Hes/Her genes inhibits the translation rate of transcribed mRNA (Murray et al. 2021). Numerical simulations were previously used to explore model behaviour and a number of experimentally testable hypotheses were defined. However, qualitative analysis of the proposed model is required in order to better characterise its behaviours and allow it to be applied in other contexts.
In this study the previous model was nondimensionalised. It was shown that the p nullcline had biologically relevant extrema if the parameter η 4 is sufficiently small. The biological interpretation of this result is that for nontrivial behaviours levels of X must be sufficiently high so as to significantly downregulate the translation rate. In order that model parameters are biologically relevant it was also found that η 2 < η 4 , i.e. the effective IC50 for transcriptional repression is an order of magnitude smaller than the IC50 for the switch of translation from low to high rates.
After performing a steady state analysis, it was shown that the model possesses a unique steady state for biologically relevant parameter values. Linear stability analysis demonstrated that the unique steady state was linearly stable when the intersection of the nullclines occurs on either the left-or right-most branches of the p nullcline. In the case where the steady state arises on the middle branch of the p nullcline, it is conditionally stable. If the mRNA is more stable than the protein (η 3 < 1) one can always identify an interval of the parameter η 1 such that the unique steady state is unstable. Upon application of the Poincare Bendixson theorem, there is therefore always a range of η 1 that yields oscillatory solutions given (provided η 3 , η 4 and η 2 are sufficiently small).
Application of numerical continuation confirmed that there is a minimal value of the parameter η 3 below which oscillatory solutions can be found. Moreover, as η 1 increases from below there is Hopf bifurcation that is either subcritical or supercritical. Notably, a previous study of isolated zebrafish PSM cells has postulated a Stuart Landau model which has a supercritical Hopf bifurcation, behaviour that is consistent with the proposed model (Webb et al. 2016).
In the limit where the rate constants associated with mRNA (η 1 and η 3 ) are chosen to be relatively small, the model behaves like a relaxation oscillator. In this case the period is approximated by assuming that the trajectory is in quasi-equilibrium on the left-and right branches of the p nullcline. Close to the local maximum of the p nullcline dynamics are relatively slow and an extra term must be accounted for that describes the time taken for protein levels to increase sufficiently so as to upregulate the translation rate.
The dimensionless equations (2) have previously been proposed as an illustrative model that describes how coupled positive and negative feedback loops give rise to 'frustrated bistability' (Krishna et al. 2009). The analysis performed heres generalises the work of Krishna et al. (2009) by considering dependence of model behaviour on the parameter η 4 as well as deriving explicit formulae for the oscillator period and bounds for the domain of oscillatory solutions. Moreover, the derivation of the model in this study differs from that of Krishna et al. (2009); here we consider regulation of the transcriptional and translational products of a single gene whilst Krishna et al. (2009) consider a model for protein-protein interaction.
The role of Notch signalling in regulating the period of the segmentation clock oscillator appear to be species dependent. In the zebrafish embryo it has been shown that levels of Notch signalling are anticorrelated with the oscillator period. When Notch signalling is increased via overexpression of Delta ligand the period decreases (Liao et al. 2016). Moreover, when levels of Notch signalling are reduced via gamma secretase treatment the period of the segmentation clock increases (Herrgen et al. 2010). In contrast, when mouse and chicken embryos are pharmaceutically treated with compounds that increase levels of the Notch intracellular domain, the period of the segmentation clock is increased (Wiedermann et al. 2015). In the proposed model intercellular coupling is not explicitly accounted for. Rather, the parameter k 1 (and hence η 1 ) can act as a proxy for levels of Notch signalling (assuming that levels of NICD regulate the maximal transcription rate). In this study it has been shown that the oscillator period can either increase or decrease with η 1 . Thus the proposed model supports the hypothesis that species-specific differences in rate constants could explain the contrasting observation of the dependence of oscillator period on levels of Notch signalling.
It is notable that oscillatory solutions of the model are permitted only if mRNA is more stable than protein. Whilst in mouse fibroblasts mRNAs have been measured to be on average approximately fives times less stable than the protein that they encode, this is not true for approximately 10% of mRNAs (Schwanhäusser et al. 2011). Gene ontology analysis associates genes that encode relatively stable mRNAs with biological processes such as tissue morphogenesis, cell proliferation, phosphorylation and positive regulation of signal transduction (Schwanhäusser et al. 2011). Moreover, direct measurement of Hes1 mRNA and protein in mouse PSM tissue yielded half lives of 24.1 and 22.3 min, respectively (Hirata et al. 2002). Additionally, in zebrafish the half life of Her7 protein has been measured to be to 3.5 min at 24 • and it has been inferred using simulations that the mRNA half life is between 2 and 6 min (Ay et al. 2013). Together, these measurements suggest that Hes1/Her7 genes encode mRNA and proteins that have similar half lives.
The relaxation oscillator analysis yields a number of experimentally testable predictions. For example, a decrease in the parameter k 5 (i.e. more stable protein) would result in a smaller oscillatory period and an increase in the amplitude of protein oscillation relative to that of the mRNA. In contrast, a decrease in parameter k 2 (i.e. making the mRNA more stable) would result in a larger period of oscillation but with an unchanged oscillation amplitude. These predictions allow for the proposed model to be distinguished from delayed negative feedback models where the oscillator period is predicted to increase linearly with both the protein and mRNA half lives Lewis (2003).
In this study a qualitative analysis has been performed on a model of a gene regulatory network in which translation as well as transcription rates are regulated by the product of a pathway. The main finding is that oscillatory solutions are possible only when: the regulation of translation rate is sufficiently large and mRNA is sufficiently more stable than protein. The qualitative analysis allows for the previous model to be applied in different biological contexts.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A Nullclines
For small p the asymptote of the p nullcline is given by the line m = p/η 4 whilst for large p it is given by the line m = p. To identify turning points of the p nullcline m 2 ( p), Eq. (4) is differentiated with respect to p and the equation is solved. Turning points thus occur at A necessary condition for unique real turning points is that an inequality that is satisfied in the intervals η 4 < 1 9 and η 4 > 1.
For p c ∈ + , a further requirement is that Hence there are two unique real positive turning points if η 4 < 1 9 .

A.1 An Approximation for the Turning Points of the p Nullcline
Given that η 4 < 1/9 then, upon applying the binomial expansion to Eq. (20), the local maximum and minimum are approximated by respectively. Note that ( p 1 , m 1 ) is a local maximum and ( p 2 , m 2 ) a local minimum of the p nullcline.

A.2 Parameter Constraints in the Case of a Non Monotonic p Nullcline
Consider the definitions of η 2 and η 4 given in Eq.
(3). The parameter ratio k 4 k 3 can be expressed in terms of dimensionless parameters in the form Hence in the case η 4 < 1, the assumption that the parameters k 3 and k 4 are positive implies that Moreover, the constraints η 2 < η 4 < 1/9 imply that k 4 k 3 > 8.
Hence the translation rate that is nonlinearly regulated must be significantly larger than the background rate k 3 in order for the p nullcline to be non-monotonic. Furthermore, using the definition of η 2 in Eq.
Recall that the quasi steady state approximation for X yields Hence for the p nullcline to be non-monotonic, X must be sufficiently large in magnitude such that it can have a a strong effect on the net translation rate.

B.1 p * is a Monotonically Increasing Function of Á 1
Differentiating Eq. (6) with respect to η 1 yields As the numerator is positive, the derivative is positive if the condition 5 p * 4 + 3 p * 2 (1 + η 2 ) − 2 p * η 1 η 2 η 3 + η 2 > 0 holds for any p * that is a solution of Eq. (22). This inequality holds in the case of a unique bounded solution to the steady state problem.

Appendix C Linear stability analysis
Positivity of the determinant of the Jacobian matrix requires Substituting for yields, after rearrangement and simplification, Noting that the inequality is therefore satisfied ∀ p > 0. Hence the Jacobian determinant is positive definite.
p † is uniquely defined so long as η 1 /η 3 > m 1 . The outward unit normal on B C is Let D represent ( p † , 0). The outward unit normal on C'D' is n C D = [1, 0]. On C'D' Finally, on D'A' the outward unit normal is Hence A B C D A defines a confined set. Therefore, by the Poincare Bendixson theorem, when the unique steady state is linearly unstable, the solution is a stable limit cycle.
Equations (2) transform to and an oscillatory solution can be approximated using a slow-fast timescale analysis.
Consider the trajectory ABCDA with coordinates On the segment AB p is assumed to be in quasi-equilibrium, i.e. m = p(1 + p 2 ) η 4 + p 2 and the ODE d p dτ = dp dm is integrated from p A to p B . After applying separation of variables the time spent on the segment AB is The integral in Eq. (27)  Given that √ η 4 , on the segment BC m is a slow variable, approximated by m ∼ m B , and p rapidly increases until the trajectory reaches the point C on the right branch of the p nullcline. On the segment CD, p is assumed to be in quasi-steady state. As p > 1 η 2 , the m dynamics are approximated by dm dτ = −η 3 m.
Thus the time spent on the segment CD is On the segment DA, dynamics are fast. m is is approximated by Thus the period is approximately given by Considering the upper bound g = g( p A ) yields the estimate Hence the period is approximated by Finally, it is noted that g < η 1 . For simplicity g is represented by the upper bound η 1 . Considering leading order in η 4 yields Eq. (17). At the local maximum of the p nullcline dp/dτ is O( √ η 4 ). Hence for ∼ √ η 4 the fast-slow analysis will become inaccurate. The time for p to increase from the local maximum of the p nullcline ( p = √ η 4 ) to the IC50 for the translation switch is approximately T B B = 1 2 √ η 4 − 1 Thus the period is approximately given by In Fig. 9 the derived estimates for the oscillator period are compared with numerical estimates. In this case η 1 ∼ O( √ η 4 ). Note that the dependence of oscillator period on η 1 is no longer monotonic. For small η 1 the period increases as η 1 decreases as a result of the progression of the solution along AB (transcription is a limiting step). However, for larger η 1 the period increases with η 1 . This effect is due to the overshoot at the local maximum of the p nullcline.