Deterministic and Stochastic Models of Arabidopsis thaliana Flowering

Experimental studies of the flowering of Arabidopsis thaliana have shown that a large complex gene regulatory network (GRN) is responsible for its regulation. This process has been mathematically modelled with deterministic differential equations by considering the interactions between gene activators and inhibitors (Valentim et al. in PLoS ONE 10(2):e0116973, 2015; van Mourik et al. in BMC Syst Biol 4(1):1, 2010). However, due to complexity of the model, the properties of the network and the roles of the individual genes cannot be deducted from the numerical solution the published work offers. Here, we propose simplifications of the model, based on decoupling of the original GRN to motifs, described with three and two differential equations. A stable solution of the original model is sought by linearisation of the original model which contributes to further investigation of the role of the individual genes to the flowering. Furthermore, we study the role of noise by introducing and investigating two types of stochastic elements into the model. The deterministic and stochastic nonlinear dynamic models of Arabidopsis flowering time are considered by following the deterministic delayed model introduced in Valentim et al. (2015). Steady-state regimes and stability of the deterministic original model are investigated analytically and numerically. By decoupling some concentrations, the system was reduced to emphasise the role played by the transcription factor Suppressor of Overexpression of Constants1 (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{SOC}1$$\end{document}SOC1) and the important floral meristem identity genes, Leafy (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{LFY}$$\end{document}LFY) and Apetala1 (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{AP}1$$\end{document}AP1). Two-dimensional motifs, based on the dynamics of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{LFY}$$\end{document}LFY and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{AP}1$$\end{document}AP1, are obtained from the reduced network and parameter ranges ensuring flowering are determined. Their stability analysis shows that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{LFY}$$\end{document}LFY and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{AP}1$$\end{document}AP1 are regulating each other for flowering, matching experimental findings. New sufficient conditions of mean square stability in the stochastic model are obtained using a stochastic Lyapunov approach. Our numerical simulations demonstrate that the reduced models of Arabidopsis flowering time, describing specific motifs of the GRN, can capture the essential behaviour of the full system and also introduce the conditions of flowering initiation. Additionally, they show that stochastic effects can change the behaviour of the stability region through a stability switch. This study thus contributes to a better understanding of the role of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{LFY}$$\end{document}LFY and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textit{AP}1$$\end{document}AP1 in Arabidopsis flowering.


Introduction
Arabidopsis thaliana is a small, annual flowering plant in the Brassicaceae (mustard) family which is a favourite model organism for plant biology research due mainly to its small size, simple genome and rapid life cycle. The transition from vegetative to reproductive development, which is an initiation of flower growth, is crucial for the life cycle of any angiosperm plant like Arabidopsis thaliana (Krizek and Fletcher 2005;Ó'Maoiléidigh et al. 2014;Wang et al. 2014) as flowering on time is a key factor to achieve reproductivity of these plants. Physiological and environmental conditions of the plant regulate the timing of transition for the optimal reproductive achievement, and their reactions are integrated into a complex GRN which monitors and regulates this transition (Kardailsky et al. 1999;Levy and Dean 1998;Wellmer and Riechmann 2010). Genes and their regulatory interactions are significant factors in biological systems at the molecular level since the understanding of their impact on each other's regulation is crucial to comprehend the response of gene disturbances on flowering time (Valentim et al. 2015). Recently, the dynamics of Arabidopsis flowering time regulation has been studied using a systems approach along with experimental data to understand the effect of the genes on flowering of Arabidopsis thaliana (Daly et al. 2009;Jaeger et al. 2013;Pullen et al. 2013;Valentim et al. 2015;Wang et al. 2014).
Numerous genes appear to be acting as flowering time regulators of Arabidopsis thaliana (Ryan et al. 2015), and different pathways have been constructed to reveal the flowering of this plant (Amasino 2010;Greenup et al. 2009;Kardailsky et al. 1999;Yant et al. 2009). This complex network of many interacting genes can be dynamically modelled using systems with many equations Valentim et al. 2015;van Mourik et al. 2010;Wang et al. 2014). In this study, we consider the deterministic dynamic model of delay differential equations (DDEs) describing the flowering of the Arabidopsis species proposed by Valentim et al. (2015). This model involves core set of gene-regulator interactions, while protein-protein interactions are not explicitly included. The model is based on a feedback loop, constructed with eight genes, where six of them are internal: Apetala1 (AP1), Leafy (LFY), Suppressor of Overexpression of Constants 1 (SOC1), Agamous-Like 24 (AGL24), Flowering Locus T (FT) and FD. The other two genes are considered as external inputs: Short Vegetative Phase (SVP) and Flowering Locus C (FLC).
System behaviour of the GRNs usually cannot be understood heuristically due to the complexity of interactions in organisms. We propose a different approach by simplifying the network and studying its behaviour. Stability analysis is used to study the properties of the GRN and threshold in flowering. Moreover, such analysis provides a reliability test and more insights into the behaviour of GRN's elements.
Our stability analysis produces conditions which include the biological parameters. Such parameters are difficult to determine from the experiment, and one of our aims was to provide specific ranges for individual coefficients that secure stable solutions. To overcome this issue of complexity, we reduce the differential equation system by decoupling some concentrations before simplifying the new system using network motifs that capture essential characteristics of the floral transition. Examples of reduced Arabidopsis thaliana GRNs can be seen in the study of Pullen et al. (2013), where a complex flowering time pathway included in the model of Jaeger et al. (2013) was simplified by focusing on essential flowering genes. Following these papers, we produce a subsystem of our network with three different motifs.
Indeed, it is known that the floral meristem identity genes have an important role to control the floral meristem specification while the flower development process is starting (Irish 2010;Levy and Dean 1998;Simon et al. 1996). Thus, this minimal regulatory network consists of the main floral meristem identity genes of Arabidopsis thaliana: AP1, LFY, FT and FD where AP1 is the dominant regulatory concentration of floral initiation with LFY in Arabidopsis thaliana (Irish 2010;Wellmer and Riechmann 2010) and has a key role between floral induction to flower formation, being a junction of flowering in the GRN ). On the other hand, FT induces flowering of Arabidopsis as an inhibitor and acts similarly with LFY. Additionally, activation tagging isolates it (Kardailsky et al. 1999). Moreover, FT and transcription factor FD affect each other in the meristem as a combined activator (Wang et al. 2014). The aim of this subsystem is to construct parameter-dependent stability conditions that reflect essential behaviour of the complex network.
Another aim of this study is to investigate the properties of the simplified Arabidopsis flowering model modified with stochastic perturbations. The motifs are reflecting the essential behaviour of the complex network and can capture the significant behaviour of the full Arabidopsis flowering model and can investigate necessary conditions (threshold values of the concentrations) for the flowering initiation. The advantage of this approach is based on the realistic description of gene effects and their interactions on flowering of Arabidopsis. New sufficient conditions of mean square stability are obtained analytically for this simplified model using Lyapunov function. Analytical and numerical investigations of the stability are performed with respect to concentrations and noise terms. This paper is organised as follows: in Sect. 2, the main features of the deterministic dynamic model of Arabidopsis flowering introduced in Valentim et al. (2015) are recalled, and analytical and numerical investigations of its steady state are both conducted. Section 3 provides a simplified deterministic model by decoupling some concentrations in the full model. A comparative numerical investigation of both models is also given. Deterministic motifs of the simplified model are presented in Sect. 4 along with an analytical investigation of their steady state and their stability. Stochastic perturbations of the motifs are investigated in Sect. 5 using Lyapunov functions to obtain sufficient conditions for their mean square stability. Finally, our concluding remarks are given in Sect. 6, while further technical information can be seen in "Appendix".

Deterministic Model
The deterministic model proposed in Valentim et al. (2015) is represented schematically in Fig. 1. Here, the transcription of FT is controlled by SV P and F LC in the leaves as shown in Fig. 1. After FT is created in the leaves, it transfers to the meristem to interact with FD. They activate the SOC1 expression together and AP1 individually (Valentim et al. 2015;Wang et al. 2014). SOC1 is activated by F T /F D, AG L24 and itself. Moreover, the expression of SOC1 is repressed by SV P and F LC in the meristem. LFY is assumed to move through a positive feedback loop with the dimerisation of AG L24 and SOC1. LFY is also a positive regulator of FD and AP1. The flowering process is determined by a direct positive input interaction among LFY and AP1. When the AP1 expression is started, the transcription variable AP1 arranges the floral transition by identifying the status of floral meristem and regulating the gene expressions comprised in flower progress (Valentim et al. 2015;Kaufmann et al. 2010). Following Valentim et al. (2015, protein and RNA levels are assumed to be linearly correlated with each other. SV P l and F LC l represent the gene expression of SV P and F LC in the leaves and SV P m and F LC m in the meristem. These four components, SV P l , F LC l , SV P m and F LC m , are independent input variables for the system which are linearly interpolated from the experimental data. These considerations led to the following system of six differential equations with one delay (Valentim et al. 2015), where the functions are defined as ( j = 1, . . . , 16 and i = 1, . . . , 10).
In system (1), the variables x i are protein concentrations, which depend on time t, and represent the genes as follows: The delayed time τ = t − appears in the equations for x 1 and x 3 . The reason for this is that FT occurs in the leaves and then moves to the meristem with some time delay , which is assumed to take less than 24 h (Valentim et al. 2015). The Hill functions γ j and κ j represent activations inhibition kinetics, respectively. The coefficient n of the Hill function γ 1 represents the cooperativity in the regulation of AP1 by LFY and is assumed to be a positive integer. The meaning of the other coefficients is provided in Table 1. Their values, estimated from experimental data using polynomial data fitting in Valentim et al. (2015), are given in Table 2. The behaviour of system (1) is simulated in Fig. 2 using the parameters in Table  2 and the experimental data used in Valentim et al. (2015). The initial conditions are taken from the experimental data. The sharp rise in AP1 from 13 to 17 days after germination can be interpreted as a predictor of flowering.
As is known from laboratory experiments (see Krämer 2015), Arabidopsis thaliana is an annual plant and its flowering is limited to approximately two to four weeks after germination. In this context, a non-trivial stable steady state can be seen as an attracting  point for the flowering process. Hence, in the next section, we turn to the analysis of the steady state of the flowering model to determine its behaviour, give conditions on its initiation and investigate the terminal stages of the flowering process.

Steady State and Stability Analysis of the Deterministic Model
Steady states of the system represent equilibrium points about which the dynamics can be studied using linear stability analysis. It helps to describe the behaviour of a delayed Table 3 Unique positive steady state for concentrations (in nm), obtained by using the parameters in Table  2 and initial values in system solution by considering the trajectories in a phase space of all dependent variables. As mentioned previously, we interpret a stable steady state as an attractor for the flowering process. Therefore, if the Arabidopsis flowering is successful, then there exists at least one strictly positive stable steady state. Hence, for DDEs of form (1), which can be presented as the equilibrium pointsx = (x 1 ,x 2 , . . . ,x 5 ,x 6 ) can be found by considering the equations f i (x 1 ,x 2 , . . . ,x 5 ,x 6 ,x 6 ) = 0, atx 6 (τ ) =x 6 . In our further consideration, we assume that the independent input variables x 7 , . . . , x 10 in system (1) are constant and equal to their initial values as given in Table 4 in "Appendix", to derive the steady states. This results intō x 2 = β 4 γ 4 (x 1 ) + β 5 γ 5 (x 3 ) + β 6 γ 6 (x 5 ) x 6 = β 12 κ 15 (x 9 )κ 16 (x 10 ) where u is a constant. Here we focus on the case n = 3, which is the value obtained by fitting experimental data in Valentim et al. (2015). It is easily seen that no trivial steady state is present whenever the constant inputs x 9 and x 10 are assumed to be nonzero. To find all equilibrium points using the assumption above, we follow the steps given in "Appendix" A.1. Eliminatingx 2 from (34) and (35), we obtain a 17thdegree polynomial equation forx 3 which we do not reproduce here. Hence, it is seen thatx 6 is obtained directly from the input concentrations, whilex 1 ,x 2 ,x 3 ,x 4 andx 5 are nonlinearly linked with each other. Using values for estimated parameters (Table 2) and the independent input variables ( Table 4 in "Appendix"), it can be seen numerically that there exists a unique positive steady state, as given in Table 3. Numerical simulation of system (1), with MATLAB R2015b, showing convergence to the steady state, is presented in Fig. 3. The time for which AP1 sees a sharp rise is in agreement with the time at which the most dramatic part of the flowering takes place, and the time for which AP1 reaches its steady state is in agreement with the ending of flowering process, which has been observed between two to four weeks in laboratory experiments (Krämer 2015; Sanda et al. 1997;Valentim et al. 2015). Our simulations show that the main features of the system behaviour would not change for different values of the input variables, apart from a slight variation in the numerical values of the steady-state concentrations.
Linearisation of the nonlinear system (1) is required to analyse the local stability of this dynamic model at its steady states (x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 ). Stability analysis is used to establish threshold conditions on the model parameters for the flowering of the plant. Therefore, we analyse the linear stability of the model in detail, and explicit conditions for local stability are formulated using the Routh-Hurwitz criterion. This gives the following theorem, for which further details can be seen in "Appendix A.2".

Theorem 1 A steady state of the nonlinear system (1) is locally asymptotically stable iff all the roots of the polynomial
have negative real parts, that is iff the following conditions are satisfied, a i > 0, i = 1, . . . , 5, a 1 a 2 a 3 + a 1 a 5 > a 2 3 + a 2 1 a 4 , and and the quantities A, B, C, . . . are defined in "Appendix A.2". Otherwise, the steady state of the system is unstable.
In summary, the conditions in Theorem 1 show that the local stability of system (1) at the steady state depends on values of parameters and concentrations. Given the high dimensionality of the parameter space, it is a difficult task to fully describe regions where stability holds. Nonetheless, it is worth noting that the delay τ does not influence stability in this particular system. No bifurcation has been numerically detected in the parameter ranges considered in this work.
To reduce the number of parameters, we now introduce a simpler system which reproduces the essential behaviour of system (1). Therefore, we performed local parameter sensitivity analysis to figure out the most important parameters in GRN (see "Appendix C"), which are β 1 , β 4 , β 5 , K 1 , K 4 , K 5 and d 1 , and all belong to the first two equations. For this purpose, we consider subsystems and analyse their stability to understand the behaviour of system (1).

Deterministic Model of the Simplified Network
The complex large regulatory network represented in (1) can be simplified while still saving its core structure. By decoupling some concentrations, it is possible to reduce the number of differential equations of the large system. One can see from the analysis in the previous section that the main contribution to the dynamics is from protein concentrations related to AP1, LFY and SOC1. Indeed, from the structure of system (1), it is seen that x 4 , x 5 and x 6 can be computed explicitly from the knowledge of x 2 , x 3 and the external outputs. Hence, we focus the analysis on these genes to investigate how they contribute to the regulation of AP1.   (5). The meaning of symbols is the same as in Fig. 1 Hence, by consideringẋ i (t) = 0 for i = 4, 5 and 6, we obtain the following system of differential equations for the variables x 1 , x 2 and x 3 : where the parameters are defined by The constant u defined in (3f) represents the constant value of the FT concentration, while U 1 and U 2 determine the effect of FT and FT-FD combination on AP1 and SOC1, respectively. The network of system (5) is described in Fig. 4. The difference with the model in Fig. 1 is that AG L24 is not involved with the external input variables SV P and F LC as they are decoupled. This network consists of three internal state variables representing SOC1, LFY and AP1, which determine the main dynamics of system (5), and two external input variables FT and FT-FD combination which are constant.
The numerical solution of the non-decoupled variables SOC1, LFY and AP1 in system (5) is compared with the numerical solution of system (1) in Fig. 5. The convergence of x 1 is affected by the constant values used for x 4 , x 5 , x 6 . Using the steady-state values, which represent the highest concentrations for these variables, leads to a slightly faster converging graph for AP1 that can be seen in Fig. 5 on the right. This result shows that decoupling some concentrations on the system can still capture the essential behaviour of the complex network for these non-constant variables.
Linear stability of the simplified model and explicit conditions for local stability at its steady states (x 1 ,x 2 ,x 3 ) are formulated using the Routh-Hurwitz criterion in Theorem 2, and further details can be seen in "Appendix A.3".
Theorem 2 A steady state of the nonlinear simplified system (5) is locally asymptotically stable iff all the roots of the polynomial have negative real parts, that is iff the conditions, a i > 0, i = 1, . . . , 3 and a 1 a 2 > a 3 , are satisfied, where

Deterministic Models of Motifs
To further reduce the complexity of system (1), we use the approach in Pullen et al. (2013) and reduce system (5) from three to two equations to understand the essential characteristics of the floral transition by considering the two components, LFY and AP1, which constitute the minimal set for enabling the transition to floral meristem (Mandel et al. 1992). Here, we model minimal regulatory networks of core components consisting of the protein concentrations for LFY, AP1, FT and FD. We consider the simplified subsystem proposed in Figure 1(b) in Pullen et al. (2013) to establish the essential characteristics of the floral transition. From system (5), one can integrate the third equation to obtain x 3 in terms of x 1 and x 2 and the various constant inputs, which also depend on FT and FD. The simplified system is represented in Fig. 4 and results from considering constant SOC1 concentrations (ẋ 3 = 0). The reason we use these four genes is: AP1 and LFY are key floral meristem identity genes in the network of Arabidopsis flowering (Irish 2010;Wellmer and Riechmann 2010) and FT induces flowering through the activation of these two genes in a feed-forward circuit (Kardailsky et al. 1999) where FD has a significant role for FT signalling. AP1 and LFY activate each other in the integration of flowering signals where they are mutual transcriptional activators (Liljegren et al. 1999). As these concentrations are key floral meristem identity genes in the network, the subsystem is based on these two genes and take into account the importance of the network activators and inhibitors. Additionally, we incorporate the action of FT-FD as a combined activator/inhibitor, as suggested in Wang et al. (2014) and Pullen et al. (2013). Ignoring the change in SOC1 concentration in the network in Fig. 4, we can redefine the simplified network as shown in Fig. 6, where LFY and AP1 represent the main dynamics of this system, and FT-FD and FT are the external input variables.
The analysis of the subsystem in Fig. 6 allows to investigate the activation and inhibition processes and provides ranges for input parameters which lead to the existence of stable solutions. Here, the effect of FT-FD and FT is described by F 1 for LFY and F 2 for AP1 in the system. Junction symbol next to AP1 shows the multiple interactions from AP1 and LFY, and γ 1 and γ 4 (Hill) functions are as in system (1). a Subsystem 1, the action is on AP1 only, b subsystem 2, the action is on LFY only, and c subsystem 3, the action is on both AP1 and LFY. a Subsystem 1, F 1 = 1, Here, F 1 and F 2 are joint inhibiting (when F i < 1) and activating (when F i > 1) constants, i = 1, 2. The variables x 1 and x 2 represent AP1 and LFY, respectively, as defined before, and we select the parameters β 1 , β 4 , K 1 , K 4 , d 1 and d 2 to be the same as previously (Table 2). This assumption relies on the statistical importance of these parameters, as determined in the sensitivity analysis. We analyse subsystem (7) in three cases including different AP1-LFY activation pathways. The first one shows the inhibition and activation of FT effect on AP1 while F 1 = 1; the second one, FT-FD effect on LFY while F 2 = 1. The third case shows the equal inhibition or activation effect of FT-FD and FT on LFY and AP1 (F 1 = F 2 ), respectively. The three realisations of the FT-FD and FT actions are given in Fig. 7.
The aim of the first and second subsystems is to analyse the effect of input variables on AP1 and LFY, individually. The third subsystem is aimed to obtain the effect of input variables when they have an equal action on both main concentrations. The parameters in Table 2 are used to investigate the behaviour of the input variables whether they play an inhibitor or an activator role.

Steady States of Motifs
The steady states (x 1 ,x 2 ) of system (7) are found by considering the right-hand side of the equations equal to zero: Here, it is easily shown that the trivial solution (x 1 ,x 2 ) = (0, 0) is a stable steady state of system (7). Although gene concentrations cannot formally be zero, the trivial steady state corresponds to a state where only small quantities are present due to nonmodelled or stochastic effects. Hence, we now focus on the non-trivial positive steady states, which can be obtained through the following process. Eliminatingx 1 from the first two equations in (8), we have where β 4 F 1 − d 2x2 > 0 as we only consider positive concentrations. This gives an upper bound for existence ofx 2 for all parameter values, Rearranging (9) implies that any non-trivial steady state forx 2 satisfies the polynomial equation where ω 1 = d 1 d 2 K 4 , ω 2 = d 2 β 1 , ω 3 = β 1 β 4 and ω 4 = d 1 d 2 K n 1 K 4 . Focusing only on positive solutions, Descartes' rule of signs indicates that (11) possesses either zero or two real positive roots, while others are complex or negative. As a consequence, system (7) has either zero or two positive steady states. We will analyse the conditions for positive real roots in further sections by using the parameter values in Table 2.

Deterministic Stability of Motifs
The dynamical subsystem (7) must have at least one stable steady state to represent the behaviour of the Arabidopsis flowering. In order to determine whether the positive equilibrium points (x 1 ,x 2 ) are locally stable, we compute the eigenvalues of the Jacobian matrix evaluated at the equilibrium points. The Jacobian matrix of systems (7) is given as, Requiring that J (x 1 ,x 2 ) is Hurwitz stable leads to the necessary and sufficient stability conditionx Combined with inequality (10), a given steady-state pointx 2 must satisfȳ in order to be stable. Full details can be seen in "Appendix (A.4)". The significance of this result is that the stability range is obtained in terms of the parameters of the system and the Hill coefficient n.

Numerical Results for Deterministic Steady States and Stability of the Motifs
Steady states are explicitly important because they offer vital knowledge on the flowering state. They can be identified by the intersection of nullclines obtained from equations (8), leading to equation (9). They are plotted in (x 1 ,x 2 ) space for the parameters in Table 2 and n = 3. The results for subsystems 1 and 2 are shown in Figs. 8 and 9. In the graphs, reference points in the plane represent the values of AP1 and LFY for specific interactions. Points where the nullclines intersect represent the steady states of the system. Non-intersecting nullclines indicate that there is no single steady state for system (7) for the given values of F 1 and F 2 .
Let us now examine the stability of steady states of LFY and AP1 for the case n = 3, for which equation (11) becomes Remembering that all coefficients ω i , F j are strictly positive, it is readily seen from Vieta's formulae that equation (14) always possesses a negative root along with either two strictly positive or complex roots. Therefore, to obtain strictly positive roots, the discriminant 3 of the cubic (14) must be positive As values of ω i , i = 1, .., 4, can be calculated from the parameters in Table 2, the discriminant only depends on the unknown values of the external input variables F i , (i = 1, 2) which represent the inhibiting (F i < 1) or activating (F i > 1) action of F T /F D. From the minimum condition of discriminant ( 3 = 0), we find the critical values of F i for the existence of such roots. The plot in the (F 1 , F 2 ) space given in Fig. 10 shows the region for the existence of positive steady states, delimited by the degeneracy condition 3 = 0 which gives rise to double roots.
Subsystem 1 (F 1 = 1). Figure 8 shows the presence of a double root at F 2 = 0.0431 from which two distinct strictly positive equilibria emanate for 0.04317 < F 2 ≤ Fig. 8 Nullclines (9) for subsystem 1 (F 1 = 1) with different values of F 2 . The red curve represents LHS of (9); the other colours represent the RHS of (9). Intersections between the red curve and other curves correspond to steady states. The thick light blue line on the right of the figure represents the numerical solution of LFY and AP1 from the original system (1) (Color figure online) max{F 2 }. Hence, when no action of FT-FD on LFY is present, the inhibition of FT on AP1 starts at the value of F 2 = 0.04317 and activation can be seen for F 2 > 1. Moreover, the behaviour of subsystem 1 is similar to system (1) for F 2 ≥ 0.04317. The best match with the numerical solution of system (1) occurs for F 2 just above 1, with a best match of the steady-state value at F 2 = 1.0476. This in turn indicates an activation action of FT on AP1.
Subsystem 2 (F 2 = 1). A similar situation is seen in this case (Fig. 9). The numerical result for this subsystem indicates that in the absence of action of FT on AP1, the inhibition of FT-FD on LFY starts at the double root F 1 = 0.05185, from which it originates one stable and one unstable positive steady states. The behaviour of subsystem 2 is similar to system (1) for F 1 ≥ 0.05185, while the best match with the numerical solution of system (1) can be seen in the activation of FT-FD on AP1 for F 1 just above 1, with a best match of the steady-state value at F 1 = 1.3445. In view of such information, we use F 1 and F 2 external input variables as an activator of the LFY and AP1 in subsystem (7) to be able to obtain a compatible behaviour with system (1).

Subsystem 3.
For subsystem 3, we distinguish two cases. In the first case, FT-FD and FT are assumed to equally inhibit/activate LFY and AP1 by using the same maximum transcription rate. A saddle-node bifurcation occurs at the value F 1 = F 2 = 0.21156; hence, two distinct positive steady states exist for all larger values, as illustrated in Fig.  Fig. 9 Nullclines (9) for subsystem 2 (F 2 = 1). The red curve represents RHS of (9); the other colours represent the LHS of (9). Intersections between red curve and other curves correspond to steady states. The thick light blue line on the right of the figure represents the numerical solution of LFY and AP1 from the original system (1) (Color figure online)

Fig. 10
Minimum condition for the existence of positive steady states in the simplified system (7) 11. The numerical solutions confirm that the actions of FT on AP1 and FT-FD on LFY do not start any interaction for the flowering of Arabidopsis until the inhibition value of F 1 = 0.21156. This situation is represented in the left-hand side of the trajectory line in Fig. 11 where there is no steady state. (The solutions of (14) forx 2 are complex.) The right-hand side of the trajectory line on this figure shows the stable and unstable steady states, indicating that the Arabidopsis flowering is in process.
In the second case, we assume F 1 and F 2 may be different from each other. In this circumstance, the best match with system (1) is for F 1 = 1.3445 and F 2 = 1.0476. These results are obtained from LHS and RHS of (8) by using the estimated parameters from Table 2 and matching the steady-state values of x 1 and x 2 from Table 3. Comparison with the solution of system (1) is given in Fig. 12, showing that subsystem 3 captures well the behaviour of the full model (1) after FT-FD and FT start activating LFY and AP1, respectively. In this case, by considering the direction of the flow dx 1 dt , dx 2 dt in the (x 1 , x 2 ) phase plane of system (7) for the obtained value F 1 and F 2 (Fig. 13), it can be explicitly seen that the trivial and non-trivial upper steady states are stable, while the lower non-trivial one is unstable. The unstable steady state can be regarded as the threshold values of the concentrations for the flowering of Arabidopsis thaliana. As a consequence, if flowering is processing for some time which means the concentrations have already reached their threshold values for the flowering, then the values of concentrations can move away from an unstable steady state and converge to a non-trivial stable one, which shows the same flowering behaviour as in the full model.
On the other hand, the initial value of the concentrations over threshold influences the flowering time of the Arabidopsis thaliana which can be seen in detail in the following figure.
As can be seen in Fig. 14, in which the AP1 value is chosen just over its threshold (0.24nM), if the initial value of LFY is lower than 1.25nM, then no flowering is observed. This is in agreement with the findings of van Dijk and Molenaar (2017), according to which mutants with knocked-down LFY may not flower. A lower threshold of LFY ≈ 1nM was estimated for the flowering of these mutants (see van Dijk and Molenaar 2017, SI Figure 9).

Fig. 14
Influence of initial LFY concentration on the flowering time of Arabidopsis thaliana, predicted using the two-dimensional motif (7) (with F 1 = 1.3445 and F 2 = 1.0476). The initial value of AP1 is chosen just over the threshold (0.24 nM). There is no flowering seen for initial LFY concentrations below 1.25 nM. The flowering time can thus be accelerated depending on the chosen initial value of LFY above this threshold

Stochastic Models of Motifs
To obtain more realistic representations of the behaviour of biological systems, it is appropriate to work with stochastic differential equations (SDEs), which can be obtained by incorporating noise terms into deterministic models. The aim of this section is to introduce and study for the first time SDEs for the behaviour of Arabidopsis flowering.
There are several ways for obtaining a SDEs model. Manninen et al. (2015Manninen et al. ( , 2006) introduced a few different approaches in their papers for incorporating stochasticity into the deterministic models. For example, stochasticity can be incorporated into reaction rates, rate constants or into concentrations by using the chemical Langevin equation. In this study, we integrate stochasticity into reaction rates, starting from the following general form of stochastic nonlinear differential equations, and consider additive (stochasticity into rate of each variable) and multiplicative (stochasticity into each reaction rate) white noise forms for G(X (t)), following Mackey and Nechaeva (1994). Both types are described and analysed in the subsections below.

Stochastic Motifs with Additive White Noise
A general Itô formulation of a system of stochastic differential equations with additive white noise form can be written as where the stochastic component GdW is added into the rate of each variable. Here, G = diag[σ 1 , · · · , σ m ] describes a nonnegative real constant diagonal matrix with parameters σ j , { j = 1, · · · , m} and W (t) is represents an m-dimensional standard Brownian motion or Wiener process over t ∈ [0, T ]. The general solution of equations (17) can be written as: The stochastic version of motif (7) can then be written as: where σ 1 , σ 2 are real constants, and W 1 and W 2 are independent standard Wiener processes with increments dW i (t) = W i (t + Δ t ) − W i (t), i = 1, 2, and each independent random variable satisfies dW i ∼ √ ΔtN (0, 1). Hill functions f 1 and f 2 are defined as and the parameters are the same as in previous sections. For numerical implementations with additive white noise, the Euler-Maruyama method with fixed time step t is used to solve this Ito SDEs model, The deterministic model (7) has three steady states: two of them are stable with a trivial and a non-trivial solution, and one is an unstable, trapped between these two stable steady states. The behaviour of this system depends on the initial conditions of the concentrations. If their initial values are lower than the unstable steady state (subthreshold value of system (7) for flowering of Arabidopsis), then system will certainly reach the trivial solution which means values are insufficient for triggering process of Arabidopsis flowering. Therefore, flowering of the Arabidopsis will not occur. If their initial values are larger than the unstable steady state, the flowering of this seed will proceed, being attracted by the non-trivial stable steady states of the concentrations.
On the other hand, the behaviour of the stochastic model (18) is more complex and depends on the initial conditions and the amount of noise in each of the concentrations. So, it is not certain whether it reaches non-trivial (passing the sub-threshold for the flowering) or trivial (non-flowering) stable equilibria, a phenomenon known as "stochastic switching" (Ullah and Wolkenhauer 2011). We show the behaviour of stochastic model (18) with a time-varying histogram to see the change of the behaviour.
The initial values are fixed as (0.2, 1.2), which lie between unstable and trivial stable steady states for the parameter values from Table 2. The implementation has been performed 100 times with a fixed constant noise of 5% (σ i = 0.05).
As can be seen from Fig. 15, stochasticity can change the behaviour of the system. The solutions are initially concentrated around the initial values, and then, they are  (18) for AP1 (left) and LFY (right) using an initial condition of (x 1 , x 2 ) = (0.2, 1.2), just below the unstable steady state and a 5% constant white noise separated into two different realisations. At the end, they converged around either trivial or non-trivial stable solutions with a considerable proportion. This shows that successful solutions for the Arabidopsis flowering can be obtained by using stochastic equations system even if the initial values are under the threshold value.
We also consider the effect of the different σ values on the stochastic system (18). If we look at the initial values of the concentrations (x init ) around the unstable steady state within 5% range, 0.95x < x init < 1.05x, we obtain the results presented in Fig. 16.

Stochastic Motifs with Multiplicative White Noise
In contrast to the previous subsection, where the possibility of successful flowering was depending only on the amount of noise terms and initial values of the concentrations, here we assume that the amplitude of noise depends on the state of the system. More precisely, stochastic perturbations of the variables around their equilibrium values are assumed to be of white noise type and proportional to the distances of AP1(x 1 ) and L FY (x 2 ) from the steady-state valuesx 1 andx 2 . The question whether the dynamical behaviour of model (7) is influenced by stochastic effects is investigated by looking at the asymptotic stochastic stability of equilibrium points.
This leads to the following stochastic differential system of the Arabidopsis flowering where again σ i are positive constants, and W i are independent standard Wiener process components, The aim is to determine the flowering domain of the stochastic motifs with multiplicative white noise. These can be obtained by using a Lyapunov function approach,  (18) with random initial condition within 5% of the unstable solution, depending on the noise parameters σ i . Blue, green and yellow dots represent success rations of the flowering process for less than 50%, between 50 and 70%, and more than 70%, respectively (Color figure online) centred at the origin or at a non-trivial steady state of the system. This allows to obtain necessary stability conditions which depend on the noise parameters σ i .
Let us show that the trivial solution x = 0 of system (20) is locally asymptotically stable in probability. Using a stochastic stability approach from Khasminskii (2011), we derive that there exists a noise-dependent domain aroundx = 0 for which asymptotic stability holds. This domain thus corresponds to non-flowering conditions for the Arabidopsis thaliana GRN modelled by system (20).

Theorem 3 The equilibrium pointx = 0 of system (20) is locally asymptotically stable in probability if the conditions
Proof 3. Letx = 0 ∈ D ⊂ R 2 be an equilibrium point of the stochastic differential equations system (20) where D is defined as a positive neighbourhood of this point. Let us define a positive definite function V where θ is a strictly positive constant. Applying a differential operator L to V (x), which is acting on the function V as gives the following expression for system (20), wherex 1 =x 2 = 0. We consider the leading terms in some positive neighbourhood around (x 1 , x 2 ) = (0, 0) of (23) up to the second-order Taylor expansion, which are Using Young's inequality in (24), By grouping x 2 1 and x 2 2 , we find The system is locally and asymptotically stable in probability if LV (x) < 0; therefore, the following inequalities are required, In particular, this implies that Then, combining inequalities (27) and (28), we find, Since all parameters are positive, this inequality can be rearranged as which shows that for any σ 1 , σ 2 satisfying σ i < √ 2d i one can choose a suitable positive value of θ such that V is a local Lyapunov function of the system. Thus, the origin is locally asymptotically stable in probability.

Conclusion
In this paper, we considered a dynamic model of Arabidopsis flowering introduced by Valentim et al. (2015). This model is reconstructed with Hill functions to emphasise the importance of these functions and their effects on the concentrations. An analytical study of the deterministic model and its steady state for the full system was performed. The stability analysis was used to establish the conditions for initiating the transition to flowering. The steady states are calculated numerically with the estimated parameters taken from Valentim et al. (2015). The analysis results have shown that the system has only one positive stable steady state and that the time for which AP1 reaches the steady state is in agreement with the observed flowering time between 20 and 30 days. The Routh-Hurwitz criterion has been used to provide local stability conditions which characterise the existence of this stable steady state; details are given in "Appendix".
Given the complexity of the system, more precise conditions have been formulated by considering subsystems which focus on the dynamics of essential elements. According to our analysis for the full system, three genes, SOC1, LFY and AP1, have a strong effect on the flowering of Arabidopsis. The network has been simplified by decoupling. Analytical solution of the simplified system is still difficult; however, it illustrates specific pathways of inhibition and activation. By using these pathways, we reconstruct three different subsystems suggested in Jaeger et al. (2013) and Pullen et al. (2013). This allowed us to derive necessary and sufficient conditions for the existence of the positive steady states of these subsystems that represent the dynamics and cooperativity of the Arabidopsis flowering time regulation system. The most important floral identity genes, AP1 and LFY, are used to investigate the flowering where they are regulating each other, and the results are confirmed by experiments (Liljegren et al. 1999). The necessary and sufficient conditions for the local stability of the deterministic model have then been determined analytically, and the stability ranges are established with the estimated parameters and compared with the numerical solutions. The numerical results have confirmed that these subsystems can capture the essential behaviour of the full model by estimating the FT-FD inhibition/activation effects on LFY and AP1, and also they help to investigate the conditions (threshold values) for the initiation of flowering, which cannot be obtained from the full model. Moreover, stochastic motifs, which are extended from the deterministic ones by adding additive and multiplicative white noise terms, have been developed to obtain more realistic description of gene effects and their interactions on the behaviour of Arabidopsis flowering. The effects of stochasticity on the steady-state regimes have been observed. The numerical solutions show that the flowering behaviour of the system does not only depend on the initial values, state variables and parameters of the stochastic system but also the amount of noise terms, where the noise can change behaviour of the stability region from non-flowering to flowering through a stability switch even if the initial values are lower than the threshold values.
Our analyses, being in a good agreement with the experimental findings, bring further insights into the roles of LFY and AP1 and provide the opportunity to explore different pathways for flowering.

A.2 Linear Stability Analysis of the Full Model
Linearisation of the nonlinear system (1) is required to analyse the local stability of this dynamic model at its steady state points (x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 ). To linearise the system with time delay, the following equation is introduceḋ to describe the behaviour in a neighbourhood of the steady-state point, where x 6 (τ )). In this equation, J 0 and J τ are Jacobians of the system with respect to non-delayed and delayed variables, respectively, The matrix form of the linearised system (36) is given as where the following notation is used, Here, γ j (x i ), for j = 1, . . . , 16 and i = 1, . . . , 6 denote derivatives of γ j with respect to x i at the steady-state point. The determinant below is introduced to obtain the characteristic equation, where I is an identity matrix. This gives the following characteristic equation, where It is clear that λ = −d 6 < 0 is a root of this characteristic equation. Thus, we now only focus on the stability of P 2 (λ) by using the Routh-Hurwitz stability criterion (Gantmacher et al. 1960). (1) is locally asymptotically stable iff all the roots of the polynomial P 2 (λ) = λ 5 + a 1 λ 4 + a 2 λ 3 + a 3 λ 2 + a 4 λ + a 5 ,

Theorem 4 A steady state of the nonlinear system
have negative real parts, that is iff the following conditions are satisfied, a 1 a 2 a 3 + a 1 a 5 > a 2 3 + a 2 1 a 4 , and (a 1 a 4 − a 5 )(a 1 a 2 a 3 + a 1 a 5 − a 2 3 − a 2 1 a 4 ) > a 5 (a 1 a 2 − a 3 ) 2 , where Otherwise, the steady state of the system is unstable.
By using the Routh-Hurwitz scheme for the 3-th degree characteristic polynomial P(λ), λ 3 1 a 2 λ 2 a 1 a 3 λ 1 b 1 0 λ 0 c 1 0 where b 1 , c 1 and the coefficients of P(λ) are as follows: b 1 = a 1 a 2 − a 3 a 1 , c 1 = a 3 , , we obtain the stability conditions of P(λ) as follows: a i > 0, i = 1, .., 3 and a 1 a 2 > a 3 . All roots will have negative real part iff the coefficients of P(λ) and b 1 are all bigger than zero.
These conditions show the validity of Theorem 2. For the parameters given in Table  2, and the initial input values in Table 4, steady state of the nonlinear simplified system (5) satisfies all conditions of Theorem 2 where a 1 = 0.9679, a 2 = 0.0943, a 3 = 0.0013 and b 1 = a 1 a 2 − a 3 = 0.093.

A.4 Stability Analysis of Deterministic Motifs
The characteristic equation of Jacobian matrix (12) is obtained as, where F 1 and F 2 have been eliminated using Eq. (8). For asymptotic stability, we require that Reλ < 0. Therefore, the necessary and sufficient conditions for local stability are tr(J ) < 0 and det(J ) > 0. The first condition is already satisfied, while the second one gives Substitutingx 1 from Eq. (8) into inequality (45), we find x n 2 > (n − 1)d 1 K n 1 K 4 (d 1 K 4 + β 1 F 2 ) .
Eliminatingx n 2 from (11), we obtain and applying equation (47) to inequality (46), gives the necessary and sufficient condition for stability as given in inequality (13).

B Tables
See Table 4.

C Parameter Sensitivity
Parameter sensitivity analysis is an approach to determine the influence of each parameters on a mathematical model output. This analysis helps to figure out the most significant parameters, which have the most important effect on the behaviour of a model. For a given dynamical systeṁ the local parameter sensitivity matrix S is obtained through the solution of an augmented system composed of the system and the sensitivity equationṡ This gives an assessment of the time-dependent influence of each parameter on the solution x. Focusing on determining which parameters most contribute to the AP1 dynamics around the flowering period, we obtain the sensitivity graphs in Figs. 17 and 18. They show that the most dominant parameters, namely β 1 , β 4 , β 5 , K 1 , K 4 , K 5 and d 1 , belong to the first two equations.