A model-based strategy to investigate the role of microRNA regulation in cancer signalling networks

Abstract

In this paper we present and discuss a model-based approach to link miRNA translational control with cell signalling networks. MicroRNAs are small regulatory RNAs that are able to regulate the activity and the stability of specific messenger RNA and have been implicated in tumour progression due to their ability to translationally regulate critical oncogenes and tumour suppressors. In our approach, data on protein–protein interactions and miRNA regulation, obtained from bioinformatics databases, are integrated with quantitative experimental data using mathematical modelling. Predictive computational simulations and qualitative (bifurcation) analyses of those mathematical models are employed to further support the investigation of such multifactorial networks in the context of cancer progression. We illustrate our approach with the C-Myc/E2F signalling network, involved in the progression of several tumour subtypes, including colorectal cancer.

Appendix

The necessity of systems theory and bifurcation analysis for investigating cell signalling networks

The nature of time delay in biochemical system models is twofold. In some cases the time delay is related to processes that take an intrinsic discrete time to be accomplished, while in other cases it is a consequence of the modelling approach used, in which complex sequences of events provoke the emergence of an apparent time delay (Nikolov et al. 2008). Time delay combined with positive/negative feedback loops (regulatory structures in which the activation of a signalling event positive/negative regulates a signalling process upstream the signalling pathway) can sometimes destabilize the stable unique equilibrium of a biochemical network and vice versa. If the time delay is large/small enough, we find that periodic solutions (sustained oscillations in the levels of the proteins and RNAs integrating the network) can arise from equilibrium by the so-called Andronov–Hopf bifurcation (Goldbeter et al. 2001; Nikolov et al. 2008; Nikolov 2008). In line with this notion, oscillations have been experimentally detected in several signalling networks showing time delays in protein expression (Nelson et al. 2004). In addition, in many cases protein–protein interactions and other biochemical processes result strongly nonlinear and show saturation at sufficiently high concentrations.

When time delay, feedback loops and non-linearity appear combined in systems like signalling networks, a suitable approach for investigating its dynamics is the systems bifurcation analysis. Bifurcation theory studies persistence and exchange of qualitative properties of dynamical systems under continuous perturbations. From the point of view of dynamic systems theory, the Hopf bifurcation theorem (Marsden and McCracken 1976; Goldbeter 2002; Nikolov et al. 2008; Nikolov et al. 2010) together with other elements of the bifurcation theory is basic analytical tools to investigate pathological conditions in biological systems (Glass and Mackey 1988; Fall et al. 2002).

A typical case is that of systems depending continuously on a single parameter. Under small variations of the parameter, the systems may lose stability, and re-stabilize near another equilibrium or a closed orbit, or a larger attractor. The type of transition can be continuous, gentle and smooth, with the new equilibrium being far from the “trivial” one (Arino et al. 2006; Shilnikov et al. 2001; Nikolov et al. 2009a). The first case corresponds to a local bifurcation, the second one is a global bifurcation. Local bifurcations can be of two types: (i) the system leaves its equilibrium state and reaches a new equilibrium state; (ii) the system goes from an equilibrium state to an invariant subset generally composed of several equilibriums and curves connecting them, closed orbits, or tori, …etc. For example, two component mechanisms with autocatalysis easily generate oscillations and bistability. They also exhibit a rich structure of bifurcations to more complicated behaviour, for example pitchfork bifurcation, saddle-node bifurcation or Takens–Bogdanov bifurcation (Fall et al. 2002; Sensse and Eiswirth 2005). The most popular and elementary situation is the Andronov–Hopf bifurcation, characterized by the onset of a closed orbit, starting near the trivial equilibrium (from focus type), which is the phase portrait of a periodic solution with a period close to some fixed number (Andronov et al. 1966; Bautin 1984; Shilnikov et al. 2001; Nikolov et al. 2009b). In this work, we will restrict our investigation to the Andronov–Hopf bifurcation.

Additional calculations performed for the analysis of time delay dynamics

Here, we consider the system (2) when all constants of model are real positive numbers. We use Andronov–Hopf bifurcation analysis when the time delay τ is a bifurcation parameter. To obtain the characteristic equation we linearize the system near the equilibrium points. Let us consider a small perturbation about the equilibrium level defined as $p = \bar{p} + x,\;v = \bar{v} + y$. Thus, we obtain:

$$ \begin{gathered} {\frac{{{\text{d}}x}}{{{\text{d}}t}}} = - \delta x + c_{1} \ell^{ - \chi \tau } x - c_{2} \ell^{ - \chi \tau } y + c_{3} \ell^{ - 2\chi \tau } x^{2} - c_{4} \ell^{ - 2\chi \tau } xy \hfill \\ - c_{5} \ell^{ - 3\chi \tau } x^{3} - c_{6} \ell^{ - 3\chi \tau } x^{2} y - c_{7} \ell^{ - 4\chi \tau } x^{4} , \hfill \\ {\frac{{{\text{d}}y}}{{{\text{d}}t}}} = k_{2} x - \gamma y, \hfill \\ \end{gathered} $$

(11)

where

$$ \begin{gathered} c_{1} = {\frac{{2\bar{p}k_{1} }}{{\Gamma_{1} }}}\left( {1 - {\frac{{2\bar{p}^{2} + \Gamma_{2} \bar{v}}}{{\Gamma_{1} }}}} \right),\quad c_{2} = {\frac{{k_{1} \Gamma_{2} \bar{p}^{2} }}{{\Gamma_{1}^{2} }}}, \hfill \\ \quad c_{3} = {\frac{{k_{1} }}{{\Gamma_{1} }}}\left( {1 - {\frac{{6\bar{p}^{2} + \Gamma_{2} \bar{v}}}{{\Gamma_{1} }}}} \right),\quad c_{4} = {\frac{{2k_{1} \Gamma_{2} \bar{p}}}{{\Gamma_{1}^{2} }}}, \hfill \\ c_{5} = {\frac{{4k_{1} \bar{p}}}{{\Gamma_{1}^{2} }}},\quad c_{6} = {\frac{{k_{1} \Gamma_{2} }}{{\Gamma_{1}^{2} }}},\quad c_{7} = {\frac{{k_{1} }}{{\Gamma_{1}^{2} }}}. \hfill \\ \end{gathered} $$

(12)

Note, that function ${\frac{1}{{\Gamma_{1} + \left[ {p\left( {t - \tau } \right)} \right]^{2} + \Gamma_{2} v\left( {t - \tau } \right)}}}$ is written in Maclaurin series and we take only the linear term. The associated transcendental characteristic equation of (11), where χ is a complex parameter, has the following form:

$$ \chi^{2} + \left( {\gamma + \delta } \right)\chi + \delta \gamma = \left( {c_{1} \chi + c_{1} \gamma - k_{2} c_{2} } \right)\ell^{ - \chi \tau } . $$

(13)

The stability of the equilibrium states depend on the sign of the real parts of the roots of (13). If both diagonal elements of the Jacobian (stability) matrix, i.e. $\left( { - \delta + c_{1} \ell^{ - \chi \tau } } \right)$ and −γ, are always negative, then tr(J) never changes sign and an Andronov–Hopf bifurcation cannot occur. If (−δ + c ₁ℓ^−χτ) and −γ are of opposite sign, then k ₂ and −c ₂ℓ^−χτ must also be of opposite sign in order det(J) to be positive. Thus, in our case, the Jacobian matrix has the typically form that produce Andronov–Hopf bifurcation, i.e.

$$ J = \left| \begin{gathered} + \quad - \hfill \\ + \quad - \hfill \\ \end{gathered} \right|. $$

(14)

Note, that mechanisms like this one are called activator–inhibitor models.

Farther, we examine Andronov–Hopf bifurcation for the system (2), using time delay as the bifurcation parameter. We rewrite (13) in terms of its real and imaginary parts as:

$$ \left| \begin{aligned} m^{2} - n^{2} + \left( {\gamma + \delta } \right)m + \delta \gamma & = \ell^{ - m\tau } \left[ {\left( {c_{1} m + c_{1} \gamma - k_{2} c_{2} } \right){ \cos }n\tau } \right. \\ & \quad + c_{1} n{ \sin }\left. {n\tau } \right], \\ 2mn + n\left( {\gamma + \delta } \right) & = \ell^{ - m\tau } \left[ {\left( {c_{1} m + c_{1} \gamma - c_{2} k_{2} } \right){ \sin }n\tau - c_{1} n{ \cos }n\tau } \right]. \\ \end{aligned} \right. $$

(15)

To find the first bifurcation point we look for purely imaginary roots $\chi = \pm in,\quad n \in R,$ of (13), i.e. we set m = 0. This substitution reduces the above two equations in the form:

$$ \left| \begin{gathered} - n^{2} + \delta \gamma = \left( {c_{1} \gamma - c_{2} k_{2} } \right){ \cos }n\tau + c_{1} n{ \sin }n\tau , \hfill \\ n\left( {\gamma + \delta } \right) = \left( {c_{1} \gamma - c_{2} k_{2} } \right){ \sin }n\tau - c_{1} n{ \cos }n\tau , \hfill \\ \end{gathered} \right. $$

(16)

One can notice that if n is a solution of (16), then so is −n. Hence, we only investigate for positive solutions n of (16). Squaring both equations in (16) and adding, we get

$$ n^{4} + \left( {\gamma^{2} + \delta^{2} - c_{1}^{2} } \right)n^{2} - \left( {c_{1} \gamma - c_{2} k_{2} } \right)^{2} + \delta^{2} \gamma^{2} = 0. $$

(17)

Because we only consider the case when system (2) is unstable for τ = 0, therefore the roots of the corresponding characteristic equation,

$$ \chi^{2} + \left( {\gamma + \delta - c_{1} } \right)\chi - \left( {c_{1} \gamma - c_{2} k_{2} } \right) + \delta \gamma = 0, $$

(18)

must have positive real parts and from the Routh–Hurwitz conditions for a square polynomial

$$ \delta \gamma - \left( {c_{1} \gamma - c_{2} k_{2} } \right) > 0\quad {\text{and}}\quad \gamma + \delta - c_{1}\,<\,0. $$

(19)

Then the left-hand side of Eq. 18 is positive for large values of n and also positive for n = 0. On the other hand, for Eq. 17 we have:

$$ n_{ \pm } = \sqrt {{\frac{{c_{1}^{2} - \gamma^{2} - \delta^{2} \pm \sqrt \Updelta }}{2}}} . $$

(20)

For both of these last expressions (17 and 19) to be real positive valued the discriminant

$$ \Updelta = \gamma^{4} + 2\left( {c_{1}^{2} - \delta^{2} } \right)\gamma^{2} - 8c_{1} c_{2} k_{2} \gamma + c_{1}^{2} \left( {c_{1}^{2} - 2\delta^{2} } \right) + 4c_{2}^{2} k_{2}^{2} + \delta^{2} , $$

(21)

must be non-negative, and $c_{1}^{2} > \gamma^{2} + \delta^{2} \mp \sqrt \Updelta$. Hence, Eq. 17 has at least one simple root and n ² is the last positive simple root of this equation. Moreover, to apply the Hopf bifurcation theorem, according to (Khan and Greenhagh 1999), the following theorem in this situation is applies:

Theorem 1

Suppose that n _b is the last positive simple root of (17). Then,in(τ _b ) = in _b is a simple root of (13) andm(τ) + in(τ) is differentiable with respect toτ in a neighbourhood ofτ = τ _b .

To establish an Andronov–Hopf bifurcation at τ = τ _b, we need to show that the following transversality condition $\left. {{\frac{{{\text{d}}m}}{{{\text{d}}\tau }}}} \right|_{{\tau = \tau_{b} }} \ne 0$ is satisfied.

Hence, if denote

$$ H\left( {\chi ,\tau } \right) = \chi^{2} + \left( {\gamma + \delta } \right)\chi + \delta \gamma - \left( {c_{1} \chi + c_{1} \gamma - k_{2} c_{2} } \right)\ell^{ - \chi \tau } , $$

(22)

then

$$ {\frac{{{\text{d}}\chi }}{{{\text{d}}\tau }}} = {\raise0.7ex\hbox{${ - {\frac{\partial H}{\partial \tau }}}$} \!\mathord{\left/ {\vphantom {{ - {\frac{\partial H}{\partial \tau }}} {{\frac{\partial H}{\partial \chi }}}}}\right.\kern-\nulldelimiterspace} \!\lower0.7ex\hbox{${{\frac{\partial H}{\partial \chi }}}$}} = {\frac{{ - \chi \left( {c_{1} \chi + c_{1} \gamma - k_{2} c_{2} } \right)\ell^{ - \chi \tau } }}{{2\chi + \gamma + \delta - c_{1} \ell^{ - \chi \tau } + \tau \ell^{ - \chi \tau } \left( {c_{1} \chi + c_{1} \gamma - k_{2} c_{2} } \right)}}}. $$

(23)

Evaluating the real part of this equation at τ = τ _b and setting χ = in _b yield

$$ \left. {{\frac{{{\text{d}}m}}{{{\text{d}}\tau }}}} \right|_{{\tau = \tau_{\text{b}} }} = {\text{Re}}\left. {\left( {{\frac{{{\text{d}}\chi }}{{{\text{d}}\tau }}}} \right)} \right|_{{\tau = \tau_{\text{b}} }} = {\frac{{n_{\text{b}}^{2} \left[ {2n_{\text{b}}^{2} + \gamma^{2} + \delta^{2} - c_{1}^{2} - 2\left( {\gamma + \delta } \right)^{2} } \right]}}{{L^{2} + I^{2} }}}, $$

(24)

where $L = \gamma + \delta + \tau_{\text{b}} \left( { - n_{\text{b}}^{2} + \delta \gamma } \right) - c_{1} { \cos }n_{\text{b}} \tau_{\text{b}}$ and $I = 2n_{\text{b}} + c_{1} { \sin }n_{\text{b}} \tau_{\text{b}} - n_{\text{b}} \tau_{\text{b}} \left( {\gamma + \delta } \right).$

Let θ = n ²_b ; then (17) reduces to

$$ g = \theta^{2} + \left( {\gamma^{2} + \delta^{2} - c_{1}^{2} } \right)\theta - \left( {c_{1} \gamma - c_{2} k_{2} } \right)^{2} + \delta^{2} \gamma^{2} . $$

(25)

Then, for $g^{\prime } \left( \theta \right)$ we have

$$ \left. {g^{\prime } \left( \theta \right)} \right|_{{\tau = \tau_{\text{b}} }} = \left. {{\frac{{{\text{d}}g}}{{{\text{d}}\theta }}}} \right|_{{\tau = \tau_{\text{b}} }} = 2\theta + \gamma^{2} + \delta^{2} - c_{1}^{2} . $$

(26)

If n _b is the least positive simple root of (17), then $\left. {{\frac{{{\text{d}}g}}{{{\text{d}}\tau }}}} \right|_{{\theta = n_{\text{b}}^{2} }}\,>\,0$. Hence,

$$ \left. {{\frac{{{\text{d}}m}}{{{\text{d}}\tau }}}} \right|_{{\tau = \tau_{\text{b}} }} = {\text{Re}}\left. {\left( {{\frac{{{\text{d}}\chi }}{{{\text{d}}\tau }}}} \right)} \right|_{{\tau = \tau_{\text{b}} }} = {\frac{{n_{\text{b}}^{2} \left[ {g^{\prime } - 2\left( {\gamma + \delta } \right)^{2} } \right]}}{{L^{2} + I^{2} }}} \,<\, 0. $$

(27)

According to the Hopf bifurcation theorem (Marsden and McCracken 1976), we define the main analytical result of this paper in the form of the following Theorem 2:

Theorem 2

If n_bis the least positive root of (17), then an Andronov–Hopf bifurcation occurs asτ passes through τ _b if only if$g^{\prime } < 2\left( {\gamma + \delta } \right)^{2}$.

Corollary 2.1

Whenτ > τ _b , then the steady state$\bar{E}$of system (2) is locally asymptotically stable under the stated conditions.

Parameter values used in the numerical analysis.

$$ \begin{gathered} \Gamma_{1} = 0.1\;\left( {\upmu {\text{M}}^{2} } \right),\quad \Gamma_{2} = 0.056\;\left( {\upmu {\text{M}}} \right),\quad \alpha \in \left( {0,0.033} \right)\;\left( {\upmu {\text{Mh}}^{ - 1} } \right), \hfill \\ \beta = 0.01\;\left( {\upmu {\text{Mh}}^{-1} } \right),\quad \delta = 0.26\;\left( {{\text{h}}^{ - 1} } \right),\quad k_{1} = 0.4\left( {\upmu {\text{Mh}}^{ - 1} } \right), \hfill \\ k_{2} = 0.3\;\left( {{\text{h}}^{ - 1} } \right),\quad \gamma = 0.02\;\left( {{\text{h}}^{ - 1} } \right). \hfill \\ \end{gathered} $$