Effects of facet joint degeneration on stress alterations in cervical spine C5–C6: A finite element analysis

1.   Introduction

Mathematical modelling of gene-regulatory networks is a relatively new area which plays an important role in systems-biology investigations. The control and coordination of large sets of genes is intrinsic to the ability of multicellular organisms to produce specific types of cells, in the proper place and at the right time during development. The perception and integration of cellular and environmental signals are essential in controlling gene expression during development. The roles of specific proteins as gene-regulatory factors are well established; in addition, recent studies of small RNAs, particularly microRNAs (miRNAs), have generated considerable excitement. Lee et al. ^[1] first described this phenomenon in $1993$, though the term microRNA was introduced only in $2001$ ^[2]. MicroRNAs are a class of small non-coding RNAs. At a post-transcriptional level, they enter the RNA interference (RNAi) pathway to regulate the expression of protein-coding genes. They regulate gene expression by blocking translation (by ribosomes) and by triggering the degradation of mRNA in, for example, both Drosophila melanogaster and Caenorhabditis elegans (see ^[3,4,5] and references therein). MicroRNAs are, therefore, key components of an evolutionarily conserved system of RNA-based gene regulation in eukaryotes. In addition, they play crucial roles in many molecular interactions, including defence against viruses and regulation of gene expression during development, cell proliferation and apoptosis. For instance, DCL1 (DICER-LIKE 1) mRNA in plants is a miRNA target, and the defects associated with dcl1 mutants include over-proliferation of meristems (which contain pluripotent stem cells), conversion of normally determinate floral meristems into indeterminate meristems, over-proliferation of embryonic suspensor cells, delayed flower timing and leaf polarity defects ^[6]. Similarly, bantam miRNA, a microRNA identified in D. melanogaster, functions to repress apoptosis and to promote cell proliferation in the developing fly, by repressing the translation of the mRNA for Hid, a key activator of programmed cell death ^[6]. It is also believed that the alterations in miRNA expression patterns might be involved in cancer development in humans ^[7]; this includes evidence for cancer-related miRNAs that regulate cellular proliferation, death and tumorigenesis in a variety of tumors, including Burkitt's lymphoma, glioblastoma, colorectal, lung and breast cancers ^[8]. MicroRNAs can also control Nodal/activin signalling in some cases ^[9]. Recently, mir-34 microRNA has been identified as a key component of the DNA damage response both in the nematode C. elegans and in human breast cancer cell lines ^[10]. In addition, studies illustrate the role of a miRNA in well-established tumour-suppressor networks associated with p53 ^[11,12,13]. Moreover, microRNAs have emerged as a class of gene-expression regulators that have also been linked to environmental stress responses, such as low temperature, high-salinity and drought. Liu et al. ^[14] identified $14$ stress-inducible miRNAs using microarray data in which the effects of abiotic stresses were surveyed in Arabidopsis thaliana. Such findings augment the current view of miRNAs as ubiquitous regulators under stress conditions.

Our goal here is to study a specific regulatory module that involves three interacting genes in C. elegans, one of which encodes a microRNA (similar genetic circuity of course also arises in other contexts): lin-4 encodes a microRNA, lin-28 a RNA-binding protein and lin-14 a transcription factor, the network being shown schematically in Figure 1. In the following paragraphs we discuss their functions in C. elegans.

Figure 1.  Network of lin-4, lin-14 and lin-28; $\rightarrow$ denotes activation and $\dashv$ repression.

DownLoad: Full-Size Img PowerPoint

In the middle of the first larval stage (L1), expression of lin-4 is first observed ^[15]. Up-regulation of this gene results in the down-regulation of the production of LIN-14 protein, which then allows the transition from expression of L1 stage to the expression of L2 (second larval) developmental events to occur ^[16].

In an experimental study, Feinbaum and Ambros ^[15] observed that over-expression of lin-4 in the L1 stage results in precocious down-regulation of LIN-14 protein and precocious expression of larval phenotypes, and noted that most of the precocious phenotypes in lin-4 over-expressing lines were strong, consistent with a potent premature repression of LIN-14 in C. elegans. However, the L1 defects (precocious expression of L2 cell lineage patterns in the L1 stage) were relatively weak. Perhaps lin-4 is not expressed at sufficiently high levels for full repression of LIN-14 synthesis until later L1. The authors of ^[15] pose a question "Does lin-4 level function as a gradient or a switch?".

In the current study, we try to give a possible answer of the question posed by ^[15] via mathematical modelling.

In summary, several heterochronic genes, such as lin-4, lin-14 and lin-28, collaborate to control the timing of specific postembryonic developmental events in C. elegans. How these three genes interact to control a particular stage-specific event of the lateral hypodermal cell lineages is described in ^[17]. Here, we also focus on these specific genes: Our goal is to develop and analyse a mathematical model for a specific gene-regulatory network (described below) in which a particular miRNA plays crucial role. It is hoped that the results will be of more general relevance in enhancing the understanding of the types of phenomena just described.

2.   Assumptions and model formulation

In developing the mathematical model, we adopt the general assumption that biochemical reactions are reversible. Next we summarize the mathematical notation and existing experimental data, along with the specific assumptions upon which the mathematical model is built (see also ^[18]).

(A1): $[N]$ and $[L]$ denote the concentrations of the mRNAs of lin-14 and lin-28, respectively, and $[\Psi]$ denotes the concentrations of the lin-4 miRNA.

(A2): Seggerson et al. ^[19] have suggested that lin-28 and lin-14 are repressed by lin-4 during normal development by a mechanism that acts on their respective mRNAs after translation initiation. This inhibition is reversible and seems to involve interference with the growing protein chain that is being translated from the mRNA. The microRNA interactions occur at the far (3') end of the mRNA molecule, even though the process of translation is initiated at the opposite (5') end. We denote the concentrations of the inactive complexes of $[N]$ and of $[L]$ with $[\Psi]$ by $[C_N]$ and $[C_L]$, respectively. These reactions are reversible and are represented by

where $k_{i}$'s are the rate constants.

(A3): We treat lin-14 and lin-28 as positively regulating each other, although in both cases this is in reality mediated through two successive negative regulatory steps, see Figure 7 of Seggerson et al. ^[19]; it is believed that the mutual positive regulation of lin-28 and lin-14 is achieved through a two-step negative regulation involving the miRNA let-7 ^{[19,20,21,22]}.

Figure 2.  The stability of the system (3.2) around $E_2(0.22, 2.02, 0.05)$ for $\zeta = 2 < {\zeta}^{[tc]} = 3.2$ and remaining parameter values given in Table 3. The initial conditions are expressed as $N_0(0)\in (0.02, 3)$, $L_0(0) = 0.136$, $\Psi_0(0) = 0.3$.

DownLoad: Full-Size Img PowerPoint

Figure 3.  Bistability of the system (3.2) for $\zeta = 4$ (that is, for a $\zeta$ lies between $\zeta^{[tc]} = 3.2$ and $\zeta^{[sn]} = 5.1$) and remaining parameter values given in Table 3. The initial conditions are provided as $N_0(0)\in (0.02, 3)$, $L_0(0) = 0.136$, $\Psi_0(0) = 0.3$.

DownLoad: Full-Size Img PowerPoint

Figure 4.  The stability of the system (3.2) around $E_1(0, 0, 0.6)$ for $\zeta = 6 > \zeta^{[sn]} = 5.1$ and remaining parameter values given in Table 3. The initial conditions are same as that of Figures 2-3.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Bistability (for $\zeta^{[tc]} = 3.2 < \zeta < \zeta^{[sn]} = 5.1$), with the transcritical bifurcation at $\zeta = {\zeta}^{[tc]} = 3.2$, and the saddle-node at $\zeta = {\zeta}^{[sn]} = 5.1$, for the parameter values given in Table 3. The solid lines represent loci of stable equilibria and the dotted lines the unstable ones.

DownLoad: Full-Size Img PowerPoint

Figure 6.  Normal Progression: The solution of the system (3.2) for $\zeta = 2.5+ 4tanh (0.001*t)$ and remaining parameter values are given in Table 3.

DownLoad: Full-Size Img PowerPoint

Figure 7.  Abnormal Progression: The solution of the system (3.2) for $\zeta = e+ dtanh (ct)$, where the parameter values are same as that of Figure 6.

DownLoad: Full-Size Img PowerPoint

(A4): lin-4 is robustly expressed from late L1 through to peak levels at L3 and persistently into adulthood ^[15]. We take lin-4 to have expression rate $\zeta$; by increasing the value of the parameter $\zeta$ we will be able to investigate its effect from the late L1 stage through to the L2 stage. Note that the numerical value of the parameter $\zeta$ is influenced by the amount of food during the time of hatching as lin-4 is repressed by starvation and activated by nutrition. In a course of time, we would be able to find that for a normal progression from L1 stage to L2 stage, an increasing value of $\zeta$ is needed (that is, a continuous nutrition is required).

A schematic diagram of the gene-regulatory network of lin-4, lin-14 and lin-28 is given in Figure 1. We model the degradation of lin-14 by the term $\beta_1[N][\Psi]$ and that of lin-28 by the term $\beta_2[L][\Psi]$.

(A5): The processes involved in the synthesis of lin-14 and lin-28 (i.e. transcription, translation and post-translational modification) are combined together, since these processes occur on a timescale of minutes to hours. The effect of a transcription factor on the transcription rate of a gene that it upregulates is taken to be described by Hill function given in the case of regulation of $N$ by $L$, for example, by $H([L]) = \frac{[L]^{m}}{[L]^{m}+\theta_{ln}^m}$. The Hill coefficient $m$ governs the steepness of the output function (the larger $m$ the more step-like the function) and we shall limit attention to the simplest Michaelis-Menten case $m = 1$ . The parameters $\theta_{ln}$ and $\lambda_{ln}$ respectively represent the threshold on a target $N$ of the regulatory influence of $L$ and the maximum rate of production of $N$ induced by $L$. Each species spontaneously decays due to degradation, at a rate $\mu_i$ ($i = 1, 2, ..., 5$).

The resulting model then reads

Tables 1 and 2 summarise the notation.

Table 1.  Variable.

Variable name	Species		Units
$[N]$	lin-14 mRNA		$\mu$M
$[L]$	lin-28 mRNA		$\mu$M
$[\Psi]$	lin-4 miRNA		$\mu$M
$[C_N]$	Complex formed by $[N]$ and $[\Psi]$		$\mu$M
$[C_L]$	Complex formed by $[L]$ and $[\Psi]$		$\mu$M

---

 | Show Table

DownLoad: CSV

Table 2.  System parameters.

Parameter	Role	Units
$\lambda_{ln}$	Maximum production rate of $[N]$ induced by $[L]$	$\mu$M h$^{-1}$
$\lambda_{nl}$	Maximum production rate of $[L]$ induced by $[N]$	$\mu$M h$^{-1}$
$\theta_{ln}$	Threshold on $[N]$ by the influence of $[L]$	$\mu$M
$\theta_{nl}$	Threshold on $[L]$ by the influence of $[N]$	$\mu$M
$k_{i}$, $\beta_{i}$,	Biochemical reaction rates ($i=1, 2.$)	$\mu$M$^{-1}$ h$^{-1}$
$k_{-i}$	Biochemical reaction rates ($i=1, 2.$)	h$^{-1}$
$\mu_i$	Decay rates ($i=1, 2, ..., 5.$)	h$^{-1}$
$\zeta$	Constitutive expression of $[\Psi]$	$\mu$M h$^{-1}$

---

 | Show Table

DownLoad: CSV

To complete the mathematical formulation of the above model we adopt for definiteness the following initial conditions:

where $N^{(0)}$, $L^{(0)}$, $\Psi^{(0)}$, ${C_N}^{(0)}$ and ${C_L}^{(0)}$ are taken to be the steady state values that arise when lin-4 is present (we mimic the transition from L1 stage to L2 stage by introducing a non-zero $\zeta$ at $t = 0$). The system (2.2) has two physiological meaningful steady states,

where $N _E$, $L_E$, $\Psi_E$, $C_{NE}$ and $C_{LE}$ are non-zero and can be obtained from an algebraic equation of order seven (see ^[18]). The steady state represented in Eq (2.5) denotes the early L1 stage where lin-14 and lin-28 are present and lin-4 is not yet expressed highly. The steady state in Eq (2.4) represents the L2 (or precocious L2) stage where lin-14 and lin-28 are absent and lin-4 is expressed highly.

We now introduce dimensionless quantities in according to

the timescale being based on the degradation rate of $[N]$, which is expected to be slow.

This gives

where the dimensionless parameters are defined by

and the initial conditions (2.3) become

We note that a linear combination of (2.6c)–(2.6e) can be chosen to eliminate the (fast) complexification terms, namely

which plays a key role in describing the dynamics in the limit $\epsilon\, \rightarrow\, 0$. We now indicate the basic assumption on parameter sizes in the system (2.6). Formation and dissociation of a complex typically operates on a timescale of the order of second, while it takes many minutes to a few hours for gene regulation to complete (transcription and translation) ^[23]; therefore we take $\epsilon\ll 1$. In the limit $\epsilon\rightarrow 0$ the order of system (2.6) is reduced and the associated singular perturbation problem is addressed in ^[18]. Here we proceed directly to the quasi-steady limit that holds for $\tau = O(1)$, the behaviour of which is significantly more accessible to analysis.

3.   The reduced system

Taking the limit $\epsilon \rightarrow 0$ with $\tau = O(1)$ in (2.6) gives

at leading order, together with

where the nonlinear functions $f_N(N_0, L_0, \Psi _0)$, $f_L(N_0, L_0, \Psi _0)$ and $f_{\Psi}(N_0, L_0, \Psi _0)$ are defined by

The reduced system (3.2) subject to the initial data

determines the behaviour for $\tau = O(1)$ and we next identify a number of its properties by a linear stability analysis. We denote a steady state solution to (3.2) by $(N_s, L_s, \Psi_s)$, so that

The system (3.2) has two kinds of steady states, of the form $E_1\equiv(0, 0, \Psi _{e}^{0})$ and $E_2 = (N_e, L_e, \Psi _e)$, where $\Psi_{e}^{0} = \frac{\zeta}{\mu_3}$, and $N _e$, $L_e$ and $\Psi_e$ are non-zero and can be obtained from a quintic equation (see ^[18]). Applying standard linear stability methods to $E_1$ and $E_2$ leads to the following propositions (details of the analysis are given in ^[18]).

Proposition 1. (ⅰ) The equilibrium point $E_1$ is stable if and only if $\zeta > {\zeta}^{[tc]}$, where

${\zeta}^{[tc]} = \frac{\mu_{{3}}(\sqrt{{ \theta_{{{\it ln}}}}^{2}{\theta_{{{\it nl}}}}^{2}{\beta_{{1}}}^{2}{\mu _{{2}}}^{2}-2\, {\theta_{{{\it ln}}}}^{2}{\theta_{{{\it nl}}}}^{2}\beta _{{1}}\mu_{{2}}\beta_{{2}}+{\theta_{{{\it ln}}}}^{2}{\theta_{{{\it nl} }}}^{2}{\beta_{{2}}}^{2}+4\, \theta_{{{\it ln}}}\theta_{{{\it nl}}} \beta_{{1}}\beta_{{2}}\lambda_{{{\it nl}}}\lambda_{{{\it ln}}}}-(\theta_{{{\it ln}}}\theta_{{{\it nl}}}\beta_{{1} }\mu_{{2}}+\theta_{{{\it ln}}}\theta_{{{\it nl}}}\beta_{{2}})) }{2\theta_{{{\it ln}}}\theta_{{{\it nl}}}\beta_{{1}} \beta_{{2}}}$.

(ⅱ) $E_1$ becomes unstable at $\zeta = {\zeta}^{[tc]}$ via a transcritical bifurcation.

Proof. (ⅰ) The Jacobian matrix of the system (3.2) around $E_1$ is given by

where

We define $C_1 = \rm{tr}(A) \times M-\det(A)$, where $M$ is the sum of the second order principal minors of $A$.

Obviously, we have $\rm{tr} A < 0$, and we find that for $\zeta > {\zeta}^{[tc]}$, $\det A < 0$ and $C_1 < 0$. Hence the proof follows from Routh-Hurwitz criteria.

(ⅱ) We notice that at $\zeta = {\zeta}^{[tc]}$, $\det A = 0$ which means that $A$ has one zero eigenvalue. Now when $\zeta = {\zeta}^{[tc]}$ the other two eigenvalues are given by $-(\mu_{{3}}+\mu_{{2}}\mu_{{3}}+\beta_{{1}}\zeta+\beta_{{2}} \zeta)/\mu_{{3}}$ and $-\mu_{{3}}$, that is both of them are negative. Let $V$ and $W$ are the eigenvectors corresponding to zero eigenvalue of the matrix $A$ and $A^T$, the transpose of $A$, respectively. We obtain that $V = (v_{1}, v_{2}, v_{3})^T$, $W = (w_{1}, w_{2}, w_{3})^T$, where

Next we see $W^T[F_{\zeta}(E_{1}, {\zeta}^{[tc]})] = 0$, $W ^T[DF_{\zeta}(E_{1}, {\zeta}^{[tc]})V]\neq 0$ and $W^T[D^2F(E_{1}, {\zeta}^{[tc]})(V, V)]\neq 0$, where $F_{\zeta} = (\frac{\partial f_N}{\partial \zeta}, \frac{\partial f_L}{\partial \zeta}, \frac{\partial f_{\Psi}}{\partial \zeta})^T$, $f_N = f_N(N, L, \Psi)$, $f_L = f_L(N, L, \Psi)$ and $f_{\Psi} = f_{\Psi}(N, L, \Psi)$ are defined in (3.3), $F_{\zeta}(E_{1}, {\zeta}^{[tc]})$ is the value of $F_{\zeta}$ evaluated at $E_1$ for $\zeta = {\zeta}^{[tc]}$, and

where $D^2F\, \in \, R^{3\times3\times3}$ is defined in Appendix A. Therefore, by the Sotomayor theorem ^[24] the system possesses a transcritical bifurcation around $E_1$, see Figure 5.

Proposition 2. (ⅰ) For ${\zeta}^{[tc]} < \zeta < {\zeta}^{[sn]}$, the system is bistable and we denote the stable states by $E_1$ and $E_2^{I}$, and the unstable one by $E_2^{II}$ (Figures 3 and 5). $E_2^{I}$ is linearly stable if $\Omega_i < 0$, $i = 1, \, 2, \, 3$, where the $\Omega_i$ are defined in the Appendix B. (ii) $E_2$ (or specifically $E_2^{I}$) experiences a saddle-node bifurcation at $\zeta = {\zeta}^{[sn]}$, where ${\zeta}^{[sn]}$ is determined from $\det \, B({\zeta}^{[sn]})$ = 0, $B$ being the Jacobian matrix of the system (3.2) at $E_2$ (the lengthy details are again given in ^[18]).

Proof. (ⅰ) We examine the bistability of the system numerically (in section 4), that when the Eq (3.4) have only two real positive roots $E_2^{I}$ and $E_2^{II}$, then the system become bistable around $E_1$ and $E_2^{I}$ (unstable around $E_2^{II}$) for $\zeta > {\zeta}^{[tc]}$; (Figures 3 and 5). However the local stability of around $E_2^{I}$ can be shown as follows:

Let us denote the Jacobian matrix of the system (3.2) around $E_2^{I}(N_e, L_e, \Psi_e)$ by $B = (b_{ij})\in R^{3\times 3}$. Here we note that the Jacobian matrix is a full matrix. Although the Routh-Hurwitz criteria can be used, but the expressions for the conditions for local stability are even more complicated than those found for the equilibrium point $E_1$. To analyse the stability of the system (3.2) around $(N_e, L_e, \Psi_e)$, we follow a different approach and get a simpler set of conditions for local stability.

To show the first claim on the asymptotic stability of the equilibrium $E_2$, we use the method of first approximation. If $B^{[2]} = \left(\bar b_{ij} \right)$ be the second computed matrix of matrix $B$, then the matrix $B^{[2]}$ is given by,

Consider a diagonal matrix $D = diag(\Psi _e, L_e, N_e)$, then $B^{[2]}$ is similar to $Q\equiv DB^{[2]}D^{-1} = (q_{ij})$, where $Q\in {\bf{R}}^{3\times 3}$ with elements

Now it is easy to see that the matrix $B^{[2]}$ is stable if and only if $Q$ is stable. Since the diagonal elements of the matrix $Q$ are negative, by the Gershgorin's theorem it is stable if it is diagonally dominant in rows. Set

and $\vartheta ^* = \max \{g_{1}, g_{2}, g_{3}\}$. Obviously, when the conditions $\Omega_1 < 0$ and $\Omega_2 < 0$ are met, then $\vartheta ^* < 0$, implying diagonal dominance. Next we have

In view of the third condition $\Omega_3 < 0$, we see $\det B < 0$, which completes the proof.

(ⅱ) Let $V$ and $U$ be the eigenvectors corresponding to eigenvalue $0$ of the matrix $B$ and its transpose, respectively. We obtain that $V = (\varphi _1, \varphi_2, \varphi_3)^T$, $U = (u_1, u_2, u_3)^T$, where $\varphi_i$ and $u_i$ are defined in Appendix C. Next we obtain $U^T[F_{\zeta}(E_{2}, {\zeta}^{[sn]})]\neq 0$, $U^T[D^2F(E_{2}, {\zeta}^{[sn]})(V, V)]\neq 0$, where $D^2F$ is given in Appendix A; therefore, system experiences saddle-node bifurcation ^[24] around the positive interior equilibrium $E_2$ at $\zeta = {\zeta}^{[sn]}$, see Figure 5.

4.   Numerical simulations

In section 1, we report that the larval development progressions of C. elegans (whether normal or abnormal) depends on the level of expression of miRNA lin-4 ^[15]. The normal development progression means the transition from expression of the L1 stage to the expression of second larval stage (L2) where as an abnormal progression refers to a premature adoption of L2 (precocious L2) by skipping the first larval development stage (L1) ^[16].

In this section we perform numerical simulations of the reduced system (3.2) to illustrate how the current mathematical modelling reflects these two types of (i.e normal and abnormal) larval development progressions. The numerical simulations are undertaken with the set of parameter values reported in Table 3 which also convenient for verifying the properties of the system determined analytically. These computations are performed using Matlab routine ode15s and XPP-Auto for stiff ODEs. The parameter values are hypothetical; however, they are closely related to the existing literature such as ^[16,23].

Table 3.  Dimensionless parameter values chosen (after removing $\tilde{}$).

Parameter	value	Parameter	value	Parameter	value
$\theta_{ln}$	7	$\lambda_{ln}$	1	$\beta_{1}$	0.01
$\theta_{nl}$	8	$\lambda_{nl}$	1	$\beta_{2}$	$0.015$
$k_{2}$	$0.8$			$\mu_2$	$0.0129$
$\mu_3$	$10$	$\mu_4$	$100$	$\mu_5$	$8$

---

 | Show Table

DownLoad: CSV

4.1.   Dynamics of the system (3.2) for constant $\zeta$

First of all we illustrate the stability of the system (3.2) around $E_2$ for $\zeta < {\zeta}^{[tc]}$, see Figure 2. Note that the stability of $E_2$ represents the early L1 larval stage where lin-4 is expressed very low and both the lin-14 and lin-28 are highly expressed. Next, we illustrate the bistability of the system (3.2) for $\zeta^{[tc]} < \zeta < \zeta^{[sn]}$: see Figure 3. Figure 3 illustrates how the system is settling to two stable states over time depending on their initial conditions. Then, we demonstrate the stability of the system (3.2) around $E_1$ for $\zeta > {\zeta}^{[tc]}$: see Figure 4. The stability of $E_1$ indicates the L2 larval stage where lin-4 is expressed very high and both the lin-14 and lin-28 are absent. Figures 2-4 are obtained by varying the values of the parameter $\zeta$ only (the other parameter values are kept same and are listed in Table 3).

Figure 5 demonstrates the transcritical bifurcation of the system (3.2) at $E_1$ for $\zeta = {\zeta}^{[tc]} = 3.2$ and the saddle-node bifurcation at $E_2$ for $\zeta = {\zeta}^{[sn]} = 5.1$.

In interpreting Figures 3 and 5, we note that for the set of parameter values listed in Table 3 and $3.2 < \zeta = 4 < 5.1$, the system (3.2) attains the stable state $E_2^{I}(0.18, 1.55, 0.09$) if the initial conditions are $N_0(0)\geq 0.023$, $L_0(0) = 0.136$ and $\Psi_0(0) = 0.3$; and it attains the stable state $E_1(0, 0, 0.4$) if $N_0(0) < 0.023$ (where $L_0(0) = 0.136$ and $\Psi_0(0) = 0.3$). Therefore, certain levels of expressions of lin-14 and lin-28 are required for L1-specific fates to occur. This means certain levels of expression of lin-14 and lin-28 are required for normal progression from the L1 stage (that is from $E_2$) through to the L2 stage (that is to $E_1$); on the other hand, if the initial levels of expression of lin-14 and lin-28 are below the required level (that is, $N_0(0) < 0.023$ and $L_0(0) < 0.136$), then the development stage L1 is skipped (as $E_2$ is not attained) resulting in a premature adoption of L2 (precocious L2) which illustrates an abnormal progression.

4.2.   Dynamics of the system (3.2) for time-dependent $\zeta$

In section 4.2, we noted that how the dynamics of the system (3.2) behave for the different constant values of $\zeta$ (constant over time $t$). It is observed that the system shows all possible scenarios (abnormal and normal progression of the larval stages) while we increase the value of $\zeta$ and (as implied by the discussion in section 1) we take an increasing value of $\zeta(t)$ to describe the up-regulation of miRNA lin-4 which causes the down-regulation of LIN-14 and LIN-28, resulting in a transition from the first developmental stage (L1) to the second developmental stage (L2) ^[16]. For these purposes we reinstate the initial data (2.4) and (2.5).

For definitions, we adopt the specific form $\zeta = A+B tanh (\sigma t)$ to describe the transition between two constant values (in keeping with its biological interpretation) with $\sigma$ sufficiently small that the system evolves in a quasi-steady fashion.

The following numerical simulations illustrate the corresponding changes in the solutions of the non-autonomous ODE system (3.2) compared to that of autonomous one. Note that the steady states $E_1$ and $E_2$ obtained earlier for a constant $\zeta$ will now be (slowly varying) functions of $t$. Figure 6 illustrates how the solution of the non-autonomous system (3.2) converges from $E_2$ (that is from L1 stage) to $E_1$ (that is to the L2 stage) over time for the initial conditions $N_0(0) = 0.25$, $L_0(0) = 2.124$, $\Psi_0(0) = 2.5$, rapidly moving between branches at $t = 4100.$ These initial conditions are taken in light of the initial conditions represented by Eq (2.5) of the original system (2.2).

On the other hand if the initial conditions are chosen in such as way which correspond the Eq (2.4) of the original system (2.2) (for example, $N_0(0) = 0$, $L_0(0) = 0$, $\Psi_0(0) = 2.5$), then $E_2$ converges (or $E_2$ is not attained) to $E_1$ very quickly, as might be anticipated. This represents an abnormal progression of the larval stages (precocious L2 to L2), see Figure 7.

Figure 8 demonstrates both the normal and abnormal progression of the larval stages in response to the time dependent $\zeta (t)$ for the both set of initial conditions used in Figures 6–7. In interpreting Figures 6-8, we again note that certain levels of expressions of lin-14 and lin-28 are necessary for a normal progression from the L1 stage to the L2 stage; otherwise, the development stage L1 is skipped and an abnormal progression takes place. These results are consistent with the results of experimental study of ^[25], where the authors studied the development of the larval stages of C. elegans based on genetic epistasis and expression analysis of heterochronic genes and observed that a high level of lin-14 is necessary for L1-specific fates to occur, while L2-specific fates occur only when both the lin-14 and lin-28 are at very low levels or off (Figure $5$ of ^[25]).

Figure 8.  This figure combines both the normal and abnormal progression of the larval stages in response to the time dependent $\zeta$ for the same set of parameter values used in Figures 6-7. One of the double curves represents the normal progression (i.e L1 to L2) and other an abnormal one (i.e Precocious L2 to L2). The current figure is the corresponding representation of Figure 5 for the case of time dependent $\zeta$.

DownLoad: Full-Size Img PowerPoint

5.   Discussion

We have proposed a mathematical model of a microRNA regulated gene network in C. elegans. The network consists two mutually activating heterochronic genes lin-14 and lin-28 and a microRNA lin-4 of which negatively regulates both of them. The reduced model (3.2) captures the qualitative behaviours, such as stability and bifurcation, of the full system (2.6). We summarize our numerical investigations of the steady state solutions in Table 4, the parameter values other than $\zeta$ being as in Table 3.

Table 4.  The table shows different sets of parameter values for which the system (3.2) has monostability and bistability.

Parameter values	Non-negative Equilibria for specific $\zeta$	Stability
$\zeta<\zeta^{[tc]} = 3.2$	Two
	For $\zeta = 2:$ $E_1$($N_e = L_e = 0$, $\Psi _e = 0.2$)	Unstable
	$E_2^{I}$ ($N_e = 0.22$, $L_e = 2.02$, $\Psi _e = 0.05$)	Stable
		(Figure 2)
$\zeta^{[tc]} = 3.2<\zeta<\zeta^{[sn]} = 5.1$	Three
	For $\zeta = 4:$ $E_1$ ($N_e = L_e = 0$, $\Psi _e = 0.4$)	Stable
	$E_2^{I}$ ($N_e = 0.18$, $L_e = 1.55$, $\Psi _e = 0.09$)	Stable
	$E_2^{II}$ ($N_e = 0.03$, $L_e = 0.25$, $\Psi _e = 0.295$)	Unstable
		(Figure 3)
$\zeta>\zeta^{[sn]} = 5.1$	One
	For $\zeta = 6:$ $E_1$($N_e = L_e = 0$, $\Psi _e = 0.6$)	Stable
		(Figure 4)

---

 | Show Table

DownLoad: CSV

It has become increasingly clear that bistability is an important recurring feature in many gene-regulatory networks. Bistability may be of particular relevance to biological systems that need switch between states. Our focus lies on two mutually activating molecular species that are down-regulated by a third species.

Interestingly, in the system (3.2), we see that bistability may arise without any cooperativity of binding, in contrast to the situation familiar for two mutually repressing genes (see ^[26] for example).

An investigation of ^[27] showed that the products of the flh-1 and flh-2 genes (encoding FLYWCH Zn finger transcription factors) function redundantly in C. elegans to repress embryonic expression of the lin-4 microRNA gene, which is normally expressed only post-embryonically. However, double mutation of flh-1 and flh-2 allows derepression of target miRNAs genes (e.g. lin-4) in embryos. Mutating the flh-3 gene (encoding a third transcription factors with a FLYWCH motif) also increases precocious expression of target miRNAs. Our model provides a mathematical interpretation for this: in the early L1 larval stage lin-4 is low, that is $\zeta$ has a sufficiently small value that can be characterised (namely, $\zeta < {\zeta}^{[tc]}$), both lin-14 and lin-28 are expressed. This corresponds to the stability of $E_2$ and the instability $E_1$ (see Figure 2 and first row of Table 4). Both lin-28 and lin-14, in turn, negatively regulate lin-29, which is believed to act as a negative regulator of early larval development, including the functions required for cell division and the expression of early larval-specific cuticle genes. Therefore, if the lin-4 microRNA expression is low, both LIN-14 and LIN-28 levels will remain high (Figure 3), and early-larval (L1-specific) stages are reiterated while later developmental events (L2-specific development) fail to occur. Conversely, if lin-4 microRNA levels remain high ($\zeta > {\zeta}^{[sn]}$), then both LIN-14 and LIN-28 levels will be low (Figure 3 and first row of Table 4), resulting in the skipping of early (L1-specific) cell lineages and developmental events, while later developmental events (L2-specific) and cell lineages occur precociously, which corresponds to stability of $E_1$ (Figure 4 and third row of Table 4). Intuitively, this means that the lin-29 gene will be turned on, resulting in the positive regulation of adult development, including the functions required for cell division and the expression for adult-specific cuticle genes; and Propositions $1$ and $2$ represent two important threshold values ${\zeta}^{[tc]}$ and ${\zeta}^{[sn]}$ of the microRNA concentration level $\zeta$. When $\zeta$ lies between these (i.e. $\zeta^{[tc]} < \zeta < \zeta^{[sn]}$), the system (3.2) has three steady states, namely $E_1$, $E_2^{I}$ (both stable) and $E_2^{II}$(unstable). Figure 5 illustrates the range of $\zeta$ for which the system (3.2) attains bistability. This switching of between two stable states takes place in between early stage L1 and the later larval stages precocious L2 which depends on the initial level of expression of lin-14, lin-28 and lin-4. This switching determines whether development progresses from L1 to L2 or whether L2-precocious patterns are reiterated (Figure 5). If the system adopts the L1 then it stays with L1 development until saddle-node bifurcation point appears (and at the saddle-node bifurcation point) whereupon the system switches to L2 development stage. This is a normal progression of the first two larval stages of C. elegans development and is consistent with the current biology ^[25]. But if the system adopts (depending on the initial level of presence of lin-14, lin-28 and lin-4) the other stable state, then this leads to the precocious adoption of later cell fates (including L2 fates appearing during L1 and so on). These events are illustrated in Figure 5, which describes an abnormal progression through these larval stages. The numerical simulations of the current model imply that a certain level of activity of the genes lin-14, lin-28 and lin-4 is necessary to allow L1-specific fates to occur (that is to adopt a normal progression), see the initial conditions of Figure 3; and this view is also supported by the experimental biology ^[25]. In summary, for a normal progression of the larval stages, three things are required (ⅰ) a certain level of activity of the genes lin-14, lin-28 and lin-4, (ⅱ) a progressive decrease in their activities owing to repression by lin-4, (ⅲ) a progressive increase in lin-4 activities (that is, increasing nutrition is required). On the other hand, the absence of lin-14 and lin-28 causes an abnormal progression through these larval stages, even though there is a progressive increase in lin-4 activities. These views are consistent with the experimental biology ^[15,16,25].

In conclusion, we have established that the simple microRNA-regulated gene network shown schematically in Figure 1 exhibits qualitative behaviour that can be interpreted in terms of key aspects of the development of C. elegans. Its partial central edge as a subnetwork of networks containing significantly more components thus warrants investigation.

Acknowledgments

The authors are grateful to Prof. Markus Owen for valuable comments. The authors are also grateful to UK-IERI for financial support to the first author.

Conflict of interest

The authors declare that no conflict of interest arises from the contents of this paper.

Appendix

A.   Notation

The notations are same as that of ^[28]. The value of $D^2F(E_{1}, {\zeta}^{[tc]})(V, V)$ is the value of $D^2F(E, \zeta)(V, V) = \left(\Upsilon_1, \Upsilon_2, \Upsilon_3\right)^T$ at $(E, \zeta) = (E_{1}, {\zeta}^{[tc]})$,

where $\Upsilon_i$'s are given by

B.   Values of $\Omega_i$

For completeness we give here expressions for the $\Omega_i$ in section 3, determined using Maple.

C.   Coefficients of the vectors $V$ and $U$.

The component of vectors $V(\varphi_1, \varphi_2, \varphi_3)$ and $U(u_1, u_2, u_3)$ are defined by