# Oscillating Shock Profiles in Relativistic Fluid Dynamics

Valentin Pellhammer  Department of Mathematics and Statistics, University of
Konstanz, 78457 Konstanz, Germany<EMAIL_ADDRESS>
###### Abstract

This note studies a model for relativistic pure-radiation fluids with
viscosity that was recently proposed by Bemfica, Disconzi and Noronha, and
shows that there are shock waves whose continuous shock profiles, if existing,
are oscillating in any variables. This behavior differs significantly from the
situation in classical fluid dynamics, in which canonical state variables are
monotone along the shock profile. This paper shows that a generic local
scenario established in [6] for a general context extends globally for this
particular model.

##  1 Introduction

As a new model for the description of viscous relativistic pure-radiation
fluids Bemfica, Disconzi and Noronha proposed in [1] a system of four partial
differential equations,

$\partial_{\beta}(T^{\alpha\beta}+\Delta T^{\alpha\beta})=0$ (1.1)

for the four velocity $u^{\gamma},\gamma=0,1,2,3$ and temperature $\theta$.
The ideal part of the energy momentum tensor is given by

$T^{\alpha\beta}=\frac{\partial(p(\theta)\psi^{\alpha})}{\partial\psi_{\beta}}=\theta^{3}p^{\prime}(\theta)\psi^{\alpha}\psi^{\beta}+p(\theta)g^{\alpha\beta},\quad
p(\theta)=\frac{1}{3}\theta^{4},$

written in Godunov variables $\psi^{\gamma}=\frac{u^{\gamma}}{\theta}$. The
dissipative part is

$\Delta
T^{\alpha\beta}=-B^{\alpha\beta\gamma\delta}(\psi)\frac{\partial\psi_{\gamma}}{\partial
x^{\delta}},$

where

$B^{\alpha\beta\gamma\delta}=\eta B_{visc}^{\alpha\beta\gamma\delta}-\mu
B_{1}^{\alpha\beta\gamma\delta}-\nu B_{2}^{\alpha\beta\gamma\delta},$

with

$B_{visc}^{\alpha\beta\gamma\delta}=\Pi^{\alpha\gamma}\Pi^{\beta\delta}+\Pi^{\alpha\delta}\Pi^{\beta\gamma}-\frac{2}{3}\Pi^{\alpha\beta}\Pi^{\gamma\delta},\quad\Pi^{\alpha\beta}=u^{\alpha}u^{\beta}+g^{\alpha\beta},$
$B_{1}^{\alpha\beta\gamma\delta}=(3u^{\alpha}u^{\beta}+\Pi^{\alpha\beta})(3u^{\gamma}u^{\delta}+\Pi^{\gamma\delta}),\quad
B_{2}^{\alpha\beta\gamma\delta}=(u^{\alpha}{\Pi^{\beta}}_{\epsilon}+u^{\beta}{\Pi^{\alpha}}_{\epsilon})(u^{\gamma}\Pi^{\delta\epsilon}+u^{\delta}\Pi^{\gamma\epsilon}).$

Additional to the viscosity parameter $\eta>0$ it depends on parameters
$\mu,\nu>0$ that were introduced to make the model causal. The model (1.1) is
a augmentation of the classical Eckart model ([7], Sec. 2.11), in order to
overcome its lack of causality. In fact, it was shown in [1, 2], that all
modes of the dissipation tensor are subluminal, if and only if

$\mu\geq\frac{4}{3}\eta,\quad\nu\leq\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1}.$
(1.2)

The non dissipative counterpart to (1.1) are given by the Euler equations

$\displaystyle\partial_{\beta}T^{\alpha\beta}=0.$ (1.3)

Fundamental solutions to the Euler equations are shock waves [3], i.e.
discontinuous solutions of the prototypical form

$\displaystyle\psi_{S}(x)=\begin{cases}\psi_{-},&x^{\beta}\xi_{\beta}<0\\\
\psi_{+},&x^{\beta}\xi_{\beta}>0\end{cases},\quad\xi^{\beta}\xi_{\beta}=1,$
(1.4)

and the question is whether these solutions correspond to solutions of the
model with dissipation (1.1). This is naturally done by means of _dissipative
shock profiles_. A shock wave (1.4) has a _dissipative shock profile_ if there
is a solution $\psi(x)=\tilde{\psi}(x^{\beta}\xi_{\beta})$ to (1.1) that
fulfills

$\displaystyle\lim\limits_{z\to\pm\infty}\tilde{\psi}(z)=\psi_{\pm}.$ (1.5)

Plugging $\tilde{\psi}$ into (1.1), integrating once and using (1.5) yields
that the shock wave (1.4) has a dissipative profile if and only if
$\tilde{\psi}$ is a solution to

$\xi_{\beta}\xi_{\delta}B^{\alpha\beta\gamma\delta}(\psi)\psi^{\prime}_{\gamma}=\xi_{\beta}T^{\alpha\beta}(\psi)-q^{\alpha},\quad
q^{\alpha}:=\xi_{\beta}T^{\alpha\beta}(\psi_{\pm})$ (1.6)

and (1.5), i.e. (1.6) admits a heteroclinic profile connecting $\psi_{-}$ with
$\psi_{+}$.

It was shown in [2] that in the case where the second inequality on (1.2)
holds strictly there are always Lax shocks that do not have a dissipative
shock profile. Based on that it was suggested that one might prefer parameters
with

$\mu\geq\frac{4}{3}\eta,\quad\nu=\left(\frac{1}{3\eta}-\frac{1}{9\mu}\right)^{-1},$
(1.7)

in the model. The choice (1.7) is referred to as _sharply causal_ , as in this
scenario there exists a mode that is luminal. The purpose of this note is to
provide further insight in the nature of shock profiles for this choice of
parameters (1.7).

In many known physical contexts shock profiles have been found to be always
monotone in natural variables. This is distinctly not the case for the present
model as the following result shows.

######  Theorem 1.

For any values $\eta,\mu,\nu>0$ satisfying (1.7), there is always a range of
Lax-shocks that have either a oscillatory profile, i.e. it is not component
wise monotone in any variables, or no profile at all.

A crucial point for the physical meaning of the shock profiles and even the
model itself is the dynamical stability. Seen from this perspective the
possible occurrence of oscillating shock profiles described in Theorem 1 is
rather delicate, as it was seen in other applications [4, 5], that oscillatory
behavior can indicate a lack of dynamical stability.

##  2 The profile equations

We follow [2] to simplify system (1.6). We assume that (1.4) is a standing
shock wave with $(1,0,0)$ as the spacial direction of propagation of the shock
front, i.e.

$(\xi_{0},\xi_{1},\xi_{2},\xi_{3})=(0,1,0,0).$

This choice causes no loss of generality by the Lorenz invariance and isotropy
of the system. With this choice (1.6) reads

$\displaystyle B^{\alpha 1\gamma 1}(\psi)\psi_{\gamma}^{\prime}=T^{\alpha
1}-q^{\alpha}.$ (2.1)

We can now consider only states $\psi$ with $\psi^{2}=\psi^{3}=0$, which
reduces (2.1) to a planar dynamical system on the state space
$\Psi=\\{\psi=(\psi^{0},\psi^{1})\in\mathbb{R}^{2}:\psi^{0}>|\psi^{1}|\\}$,
where the indices in (2.1) run over $0$ and $1$ instead from $0$ to $3$. The
temperature and the velocity are now given as

$\displaystyle\theta=(-\psi_{\gamma}\psi^{\gamma})^{-\frac{1}{2}},\quad(u,v):=(u^{0},u^{1})=\theta(\psi^{0},\psi^{1})$
(2.2)

with $u^{2}=1+v^{2}$.

######  Lemma 1.

System (2.1) has more than one rest point if and only if

$q^{1}>0,\quad(q^{1})^{2}<(q^{0})^{2}<\frac{2}{\sqrt{3}}(q^{1})^{2}.$

In this case it has exactly two rest points. These states form a standing Lax-
shock with respect to the Euler equations in right-moving or left-moving flow
if $q^{0}>0$ or $q^{0}<0$ respectively.

###### Proof.

Cf. [2], Lemma 1. ∎

As system (2.6) is invariant under the reflection $q^{0}\mapsto-q^{0}$,
$u^{1}\mapsto-u^{1}$ we can assume w.l.o.g.

$q^{0}>0,\quad q^{1}>0$

The pair of equilibria existing by Lemma 1 depend on

$\tilde{q}:=\left(\frac{q^{1}}{q^{0}}\right)^{2}$

and is given by
$(\psi_{-},\psi_{+})=(\psi_{-}(\tilde{q}),\psi_{+}(\tilde{q}))$ with

$\displaystyle\psi_{\pm}(\tilde{q})=\left(\frac{4}{3}v^{2}+\frac{1}{3}\right)^{\frac{1}{4}}\begin{pmatrix}(1+v^{2})^{\frac{1}{2}}\\\
v\end{pmatrix}\bigg{|}_{v=v_{\pm}(\tilde{q})},\quad
v_{\pm}^{2}(\tilde{q})=\frac{(2\tilde{q}-1)\mp\sqrt{\tilde{q}(4\tilde{q}-3)}}{4(1-\tilde{q})}.$
(2.3)

They exist if and only if

$\displaystyle\frac{3}{4}<\tilde{q}<1.$

For later use we note that

$\displaystyle v_{+}^{2}\in\left(\frac{1}{8},\frac{1}{2}\right).$ (2.4)

We can furthermore assume w.l.o.g. $q^{1}=1$ and hence
$\tilde{q}=(q^{0})^{-2}$, cf. [2]. Note that the limit
$\tilde{q}\to\frac{3}{4}$ refers to the zero amplitude limit of the shock wave
and $\tilde{q}\to 1$ refers to the infinite amplitude of the shock wave.

Consider parameters with (1.7). A scaling of the independent variable in (2.1)
according to $t\mapsto\frac{t}{\mu}$ and introducing the parameter

$\displaystyle\epsilon:=\frac{4\nu}{3\mu}$ (2.5)

results in the two parameter dynamical system

$\displaystyle B^{\\#}(\psi,\epsilon)\psi^{\prime}=F(\psi,\tilde{q})$ (2.6)

with

$\displaystyle F(\psi,q):=\begin{pmatrix}F^{0}\\\
F^{1}\end{pmatrix}:=\begin{pmatrix}-\frac{4}{3}\theta^{4}vu+q^{0}\\\
\theta^{4}(\frac{4}{3}v^{2}+\frac{1}{3})-q^{1}\end{pmatrix}$

and

$B^{\\#}(\psi,\epsilon):=\epsilon
B_{visc}(\psi)-B_{1}(\psi)-\frac{9\epsilon}{4-\epsilon}B_{2}(\psi),$

where

$B_{visc}(\psi)=\begin{pmatrix}u^{2}v^{2}&-u^{3}v\\\
-u^{3}v&u^{4}\end{pmatrix},\quad
B_{1}(\psi)=\begin{pmatrix}16u^{2}v^{2}&-4uv(4v^{2}+1)\\\
-4uv(4v^{2}+1)&(4v^{2}+1)^{2}\end{pmatrix},$
$B_{2}(\psi)=\begin{pmatrix}(u^{2}+v^{2})^{2}&-2(u^{2}+v^{2})uv\\\
-2(u^{2}+v^{2})uv&4u^{2}v^{2}\end{pmatrix}.$

With (2.5), system (2.6) fulfills the condition (1.7) if and only if

$\displaystyle 0<\epsilon\leq 1.$

######  Lemma 2.

Consider system (2.6) and assume $0<\epsilon\leq 1$ as well as
$\frac{3}{4}<\tilde{q}<1$. Then $\psi_{-}(\tilde{q})$ is a hyperbolic saddle
and $\psi_{+}(\tilde{q})$ is an attractor.

###### Proof.

Fix $0<\epsilon\leq 1$ and $\frac{3}{4}<\tilde{q}<1$. It was already shown in
[2], that $\psi_{-}(\tilde{q})$ is a saddle and
$\det(DF(\psi_{+}(\tilde{q})))<0$ for each choice of
$\frac{3}{4}<\tilde{q}<1$. As

$\det(B^{\sharp}(\psi,\epsilon))=\frac{9\epsilon((8+\epsilon)v^{2}+\epsilon-1)}{\epsilon-4},$

we find $\det(B^{\\#}(\psi,\epsilon))<0$ whenever
$v^{2}>\frac{1-\epsilon}{8+\epsilon}$. Looking at (2.4) and observing
$\frac{1-\epsilon}{8+\epsilon}\in[0,\frac{1}{8})$ for all $0<\epsilon\leq 1$
yields that $\det(B_{\mu}^{\sharp}(\psi_{+}(\tilde{q})))<0$ and therefore

$\det\left(B^{\sharp}(\psi_{+}(\tilde{q}),\epsilon)^{-1}\;DF(\psi_{+}(\tilde{q}))\right)>0.$

This shows that $\psi_{+}(\tilde{q})$ is an attractor or a repeller. A simple
calculation yields

$trace\left(B^{\sharp}(\psi,\epsilon)^{-1}DF(\psi)\right)=-\frac{3v}{(4-\epsilon)\det(B^{\sharp}(\psi,\epsilon))}((8+\epsilon^{2})v^{2}+\epsilon^{2}-6\epsilon-4)$

which is negative if
$0<v^{2}<\frac{-\epsilon^{2}+6\epsilon+4}{\epsilon^{2}+8}$. From the fact that
$\frac{-\epsilon^{2}+6\epsilon+4}{\epsilon^{2}+8}\in(\frac{1}{2},1]$ for all
$0<\epsilon\leq 1$ and (2.4) we can deduce

$trace\left(B^{\sharp}(\psi_{+}(\tilde{q}),\epsilon)^{-1}DF(\psi_{+}(\tilde{q}))\right)<0,$

which shows that $\psi_{+}(\tilde{q})$ is an attractor. ∎

We study (2.6) as a dynamical system dependent on two parameters in the
parameter space

$\Omega:=\left\\{(\epsilon,\tilde{q}):0<\epsilon\leq
1,\frac{3}{4}<\tilde{q}<1\right\\},$

since it covers all possible shocks and choices of parameters with (1.7) for
the model.

##  3 Rest points with non real eigenvalues

Lemma 2 makes sure that $\psi_{+}$ is an attractor for each choice of
$\frac{3}{4}<\tilde{q}<1$. As such it is either a stable focus or a stable
node, depending on weather the eigenvalues of the lineraziation of (2.1) have
a nonzero imaginary part or not. A _scaled_ version of the linarization around
a state $\bar{\psi}$ is given by

$\displaystyle B^{\\#}(\bar{\psi},\epsilon)\psi^{\prime}=A(\bar{\psi})\psi$
(3.1)

with

$A(\psi):=\begin{pmatrix}v(6v^{2}+5)&-u(6v^{2}+1)\\\
-u(6v^{2}+1)&3v(2v^{2}+1)\end{pmatrix}.$

Depending on the two parameters $(\epsilon,\tilde{q})\in\Omega$ the
eigenvalues of the matrix
$B^{\\#}(\psi_{+}(\tilde{q}),\epsilon)^{-1}A(\psi_{+}(\tilde{q}))$ determine
if $\psi_{+}(\tilde{q})$ is a stable node or a stable spiraling focus. We
therefore consider the discriminant of the characteristic polynomial of the
matrix $B^{\\#}(\psi,\epsilon)^{-1}A(\psi)$, which is given by

$\displaystyle
D(\psi,\epsilon):=trace(B^{\\#}(\psi,\epsilon)^{-1}A(\psi))^{2}-4\det(B^{\\#}(\psi,\epsilon)^{-1}A(\psi))=\frac{1}{\det(B^{\\#}(\psi,\epsilon))^{2}}\tilde{P}(\psi,\epsilon)$

with

$\tilde{P}(\psi,\epsilon)=trace(adj(B^{\\#}(\psi,\epsilon))A(\psi))^{2}-4\det(B^{\\#}(\psi,\epsilon))\det(A(\psi))),$

where $adj(B^{\\#})$ denotes the adjunct matrix of $B^{\\#}$. A simple
calculation yields

$\displaystyle\det(B^{\\#}(\psi,\epsilon))=\frac{9\epsilon}{\epsilon-4}((\epsilon+8)v^{2}+\epsilon-1),\quad\det(A(\psi))=2v^{2}-1$

and

$trace(adj(B^{\\#}(\psi,\epsilon))A(\psi))=\frac{3v}{\epsilon-4}((8+\epsilon^{2})v^{2}+\epsilon^{2}-6\epsilon-4),$

from which we deduce

$\tilde{P}(\psi,\epsilon)=\frac{9}{(\epsilon-4)^{2}}P(v^{2},\epsilon)$

with the Polynomial $P$ given by

$P(z,\epsilon)=a_{3}(\epsilon)z^{3}+a_{2}(\epsilon)z^{2}+a_{1}(\epsilon)z+a_{0}(\epsilon)$

where

$\displaystyle a_{0}(\epsilon)=\epsilon(4\epsilon^{2}-20\epsilon+16),\quad
a_{1}(\epsilon)=\epsilon^{4}-16\epsilon^{3}+84\epsilon^{2}-112\epsilon+16$
$\displaystyle
a_{2}(\epsilon)=2\epsilon^{4}-20\epsilon^{3}-24\epsilon^{2}+160\epsilon-64,\quad
a_{3}(\epsilon)=\epsilon^{4}+16\epsilon^{2}+64.$

For each $\epsilon>0$ and $\psi\in\Psi$, the sign of $D(\psi,\epsilon)$
coincides with the sign of $P(v^{2},\epsilon)$ with $v$ given according to
(2.2).

######  Lemma 3.

Considers $P$ as a parameter dependent polynomial in $z$ with parameter
$0\leq\epsilon\leq 1$. For each $0\leq\epsilon\leq 1$, $P$ has three real
zeros $w_{1}(\epsilon)<w_{2}(\epsilon)\leq w_{3}(\epsilon)$ with
$w_{2}(\epsilon)=w_{3}(\epsilon)$ if and only if $\epsilon=0$ and the
following properties.

  1. (i)

$w_{1}(\epsilon)<0$ for all $0<\epsilon\leq 1$ and $w_{1}(0)=0$.

  2. (ii)

for
$\hat{\epsilon}:=\frac{2}{3}\left(3\sqrt{6}-2\sqrt{16-6\sqrt{6}}-4\right)\approx
0.7103$ one has $w_{2}((0,\hat{\epsilon}))\subset(\frac{1}{8},\frac{1}{2})$
and $w_{2}((\hat{\epsilon},1))\subset(0,\frac{1}{8})$ as well as
$w_{2}(0)=\frac{1}{2}$,$w_{2}(\hat{\epsilon})=\frac{1}{8}$ and $w_{2}(1)=0$.

  3. (iii)

$w_{3}(\epsilon)\in(\frac{1}{3},\frac{1}{2})$ for $0<\epsilon<1$ and
$w_{3}(0)=\frac{1}{2}$, $w_{3}(1)=\frac{1}{3}$.

###### Proof.

The discriminant of $P$ as a polynomial in $z$ is given by
$D(\epsilon)=1296\epsilon^{5}(4-\epsilon)^{3}\tilde{D}(\epsilon)$ where

$\tilde{D}(\epsilon)=-4\epsilon^{5}+179\epsilon^{4}-844\epsilon^{3}+880\epsilon^{2}-32\epsilon+64.$

Setting $\zeta=\frac{1}{\epsilon}$ yields

$\tilde{D}(\epsilon)=\zeta^{-5}\tilde{D}_{2}(\zeta),\quad\text{with}\quad\tilde{D}_{2}(\zeta)=64\zeta^{5}-32\zeta^{4}+880\zeta^{3}-844\zeta^{2}+179\zeta-4.$

The polynomial

$\tilde{D}_{2}(x+1)=64x^{5}+288x^{4}+1392x^{3}+2244x^{2}+1323x+243$

has only positive coefficients and can thus not have a non negative real root.
Therefore $\tilde{D}_{2}$ has no root in $[1,\infty)$ and has a constant sign
in this interval. Since the leading coefficient of $\tilde{D}_{2}$ is positive
one has $\tilde{D}_{2}(\zeta)>0$ for all $1<\zeta<\infty$ and we can conclude
that $D(\epsilon)>0$ for all $0<\epsilon<1$. This shows the existence of three
distinct real roots $w_{1}(\epsilon)<w_{2}(\epsilon)<w_{3}(\epsilon)$, which
depend continuously on $\epsilon$.

One has

$P(z,1)=27z(3z^{2}+2z-1),$

and therefore $w_{1}(1)=-1,w_{2}(1)=0$ and $w_{3}(1)=1/3$.

Since the leading coefficient of $P$ is positive and
$P(0,\epsilon)=4\epsilon(\epsilon^{2}-5\epsilon+4)>0$ for all $0<\epsilon\leq
1$ it follows that $P$ has a negative real root, i.e. $w_{1}(\epsilon)<0$ for
all $0<\epsilon\leq 1$. This shows assertion (i).

As $P$ is cubic in $z$ with a positive leading coefficient and
$P(\frac{1}{2},\epsilon)=\frac{9}{8}\epsilon^{2}(\epsilon-4)^{2}>0$ we know
that $w_{3}(\epsilon)\neq\frac{1}{2}$ for all $0<\epsilon\leq 1$. Since
additionally $w_{3}(1)=\frac{1}{3}<\frac{1}{2}$ we can deduce that
$w_{3}(\epsilon)<\frac{1}{2}$ for all $0<\epsilon\leq 1$ by the continuity of
$w_{3}$.

One has
$P(\frac{1}{3},\epsilon)=\frac{16}{27}(\epsilon-1)^{2}(\epsilon^{2}-4\epsilon+1)$
and thus $P(\frac{1}{3},\epsilon)=0$ for a $0<\epsilon\leq 1$ if and only if
$\epsilon\in\\{2-\sqrt{3},1\\}$ and it follows that
$w_{2}(2-\sqrt{3})=\frac{1}{3}<w_{3}(2-\sqrt{3})$. Since additionally
$w_{3}(\epsilon)=\frac{1}{3}$ if and only if $\epsilon=1$ we follow by
continuity of $w_{3}$ that $\frac{1}{3}<w_{3}(\epsilon)$ for all
$0<\epsilon<1$ which thus shows (ii).

For $\epsilon>0$ one has

$\displaystyle P\left(\frac{1}{8},\epsilon\right)$
$\displaystyle=\frac{9}{512}(9x^{4}+96x^{3}-560x^{2}+256x+64)$
$\displaystyle=\frac{9\epsilon^{2}}{512}\left(9\epsilon^{2}+96\epsilon^{1}-560+256\frac{1}{\epsilon}+64\frac{1}{\epsilon^{2}}\right)$
$\displaystyle=\frac{9\epsilon^{2}}{512}\left(9\left(\epsilon+\frac{8}{3\epsilon}\right)^{2}+96\left(\epsilon+\frac{8}{3\epsilon}\right)-608\right)$

and hence $P(\frac{1}{8},\epsilon)=0$ can easily be solved for
$y:=(\epsilon+\frac{8}{3\epsilon})$. Doing this yields that
$P(\frac{1}{8},\epsilon)=0$ for a $0<\epsilon\leq 1$, i.e.
$w_{2}(\epsilon)=\frac{1}{8}$ if and only if
$\epsilon=\hat{\epsilon}:=\frac{2}{3}\left(3\sqrt{6}-2\sqrt{16-6\sqrt{6}}-4\right)$.
Since $w_{2}(1)=0$ the continuity of $w_{2}$ yields
$w_{2}(\epsilon)<\frac{1}{8}$ for each $\hat{\epsilon}<\epsilon\leq 1$. ∎

The following Theorem reveals the transition of $\psi_{+}$ from a stable node
to a stable focus for varying parameters in detail. In the case where
$\psi_{+}$ is a focus a possibly existing orbit connecting $\psi_{-}$ with
$\psi_{+}$ has a oscillating character due to its spiraling around $\psi_{+}$
at that end state, which implies Theorem 1.

######  Theorem 2.

The parameter range $\Omega$ decomposes into two open subsets $\Omega_{node}$
and $\Omega_{focus}$ and two seperatrices $Q_{1}$ and $Q_{2}$ such that the
equilibrium $\psi_{+}$ is a node if and only if
$(\epsilon,\tilde{q})\in\Omega_{node}$ and $\psi_{+}$ is a stable focus if and
only if $(\epsilon,\tilde{q})\in\Omega_{focus}$. The sepreatrices $Q_{1}$ and
$Q_{2}$ are the graphs of two functions $\tilde{q}_{1}$ and $\tilde{q}_{2}$
given by

$\tilde{q}_{1}:(0,1],\mapsto(\tfrac{3}{4},1),\quad\tilde{q}_{2}:(0,\hat{\epsilon}),\mapsto(\tfrac{3}{4},1).$

$\Omega_{focus}$ is the open subset contained between the two curves and
$\Omega_{node}$ consists of the two open sets $\Omega_{node}^{1}$ below
$Q_{1}$ and $\Omega_{node}^{2}$ above $Q_{2}$, cf. Figure 1. Furthermore

$\displaystyle\tilde{q}_{1}(1)=\frac{49}{64},\quad\lim\limits_{\epsilon\to\hat{\epsilon}}\tilde{q}_{2}(\epsilon)=1,\quad\lim\limits_{\epsilon\to
0}\tilde{q}_{1}(\epsilon)=\lim\limits_{\epsilon\to
0}\tilde{q}_{2}(\epsilon)=\frac{3}{4}$ (3.2)

###### Proof.

Consider the discriminant of the characteristic polynomial of the
linearization at $\psi_{+}(\tilde{q})$ given by

$D(\psi_{+}(\tilde{q}),\epsilon)=\frac{9}{\det(B)^{2}(\epsilon-4)^{2}}P(v_{+}^{2}(\tilde{q}),\epsilon).$

and the function

$V_{+}:\left(\tfrac{3}{4},1\right)\to\left(\tfrac{1}{8},\tfrac{1}{2}\right),\quad\tilde{q}\mapsto
v_{+}^{2}(\tilde{q})$

that links the value $\tilde{q}$ to the value for $v^{2}_{+}(\tilde{q})$ of
the equilibrium $\psi_{+}(\tilde{q})$ via (2.3). The function $V_{+}$ is onto
and strictly decreasing and has thus a continuous inverse $V_{+}^{-1}$ that is
also strictly decreasing. We define

$\displaystyle\tilde{q}_{1}(\epsilon):=(V_{+}^{-1}\circ
w_{3})(\epsilon),\epsilon\in(0,1],\quad\tilde{q}_{2}(\epsilon):=(V_{+}^{-1}\circ
w_{2})(\epsilon),\epsilon\in(0,\hat{\epsilon})$

where $w_{2}$ and $w_{3}$ are the functions defined in Lemma 3. Since
$V_{+}^{-1}$ is monotone one has $\tilde{q}_{1}<\tilde{q}_{2}$ on
$(0,\hat{\epsilon})$ and the assertions (3.2) easily follow using Lemma 3.

Since $P$ is cubic and the leading coefficient of $P$ is positive, Lemma 3
shows that $P(z,\epsilon)<0$ for $z\in(\frac{1}{8},\frac{1}{2})$ and
$\epsilon\in(0,1]$ if and only if either $\epsilon\in(0,\hat{\epsilon})$ and
$z\in(w_{2}(\tilde{q}),w_{3}(\tilde{q}))$ or $\epsilon\in[\hat{\epsilon},1)$
and $z>w_{2}(\tilde{q})$.

For $\epsilon\in(0,\hat{\epsilon})$ and
$q_{1}(\epsilon)<\tilde{q}<q_{2}(\epsilon)$ one has
$v_{+}^{2}(\tilde{q})\in(w_{2}(\tilde{q}),w_{3}(\tilde{q}))$ and therefore
$D(\psi_{+}(\tilde{q}),\epsilon)<0$ such that the eigenvalues of the
linearization (3.1) have a nonzero real part and $\psi_{+}(\tilde{q})$ is a
stable focus. If $\epsilon\in(0,\hat{\epsilon})$ and either
$\tilde{q}<q_{1}(\epsilon)$ or $q_{2}(\epsilon)<\tilde{q}$ one has
$v_{+}^{2}(\tilde{q})\notin(w_{2}(\tilde{q}),w_{3}(\tilde{q}))$ and thus
$D(\psi_{+}(\tilde{q}),\epsilon)>0$. In this case the eigenvalues of (3.1)
have a non negative imaginary part and $\psi_{+}(\tilde{q})$ is a stable node.

In the case $\epsilon\in(\hat{\epsilon},1]$ the same arguments yield that
$\psi_{+}(\tilde{q})$ is a focus if $q_{1}(\epsilon)<\tilde{q}$ and a node if
$q_{1}(\epsilon)<\tilde{q}$. ∎

Figure 1: The parameter space for the equation (2.6). $\tilde{q}$ controls the
strength of the shock and $\epsilon$ refers to the choice of dissipation. The
rest point $\psi_{+}$ is a stable focus for
$(\epsilon,\tilde{q})\in\Omega_{focus}$ and a stable node for
$(\epsilon,\tilde{q})\in\Omega_{node}=\Omega_{node}^{1}\cup\Omega_{node}^{2}$.
For every choice of sharply causal dissipation there are shock waves, whose
profiles, if they exist, can only be oscillatory.

## References

  * [1] F. S. Bemfica, M. M. Disconzi, J. Noronha: Causality and existence of solutions of relativistic viscous fluid dynamics with gravity, Phys. Rev. D 98 (2018), 104064.
  * [2] H. Freistühler: Nonexistence and existence of shock profiles in the Bemfica-Disconzi-Noronha model, Phys. Rev. D 103 (2021), 124045.
  * [3] P. D. Lax: Hyperbolic systems of conservation laws II, Comm. Pure Appl. Math. 10 (1957), 537–566.
  * [4] T. Li, H. Liu, L. Wang: Oscillatory traveling wave solutions to an attractive chemotaxis system, J. Differ. Eqs. 261 (2016), 7080 –7098.
  * [5] R. L. Pego, P. Smereka, M. I. Weinstein: Oscillatory instability of traveling waves for a KdV-Burgers equation Physica D, 67 (1993), 45 –65.
  * [6] V. Pellhammer: A generically singular type of saddle node bifurcation that occurs for relativistic shock waves, in preparation.
  * [7] S. Weinberg: Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity, Wiley, New York, 1972.


# Flexible and robust networks

S. Vakulenko1 and O. Radulescu2
1 Saint Petersburg State University of Technology and Design, St.Petersburg,
Russia,
2 DIMNP UMR CNRS 5235, University of Montpellier 2, Montpellier, France.

Abstract  We consider networks with two types of nodes. The $v$-nodes, called
centers, are hyperconnected and interact one to another via many $u$-nodes,
called satellites. This centralized architecture, widespread in gene networks,
possesses two fundamental properties. Namely, this organization creates
feedback loops that are capable to generate practically any prescribed
patterning dynamics, chaotic or periodic, or having a number of equilibrium
states. Moreover, this organization is robust with respect to random
perturbations of the system.

##  1 Introduction

Flexibility and robustness are important properties of biological systems.
Flexibility means the capacity to adapt with respect to changes of environment
whereas robustness is the capacity to support homeostasis in spite of
environmental changes. Intriguingly, it seems that biological systems could be
in the same time robust and flexible. Development of an organism is robust to
variations of initial conditions and environment, species can diversify in
order to better satisfy constraints imposed by a varying environment.

We discuss here flexibility and robustness problems for genetical networks of
a special topological structure as a model for flexible and robust systems. In
these networks, highly connected hubs play the role of organizing centers
(centralized networks). The hubs receive and dispatch interactions. Each
center interacts with many weakly connected nodes (satellites). Similar ideas,
that such a ”bow-tie” connectivity can play a role in robustness, have been
proposed by (Zhao et al. 2006). In the field of random boolean networks
(Kauffman 1969), the phase transitions from chaotic to frozen (robust) phases
were related to scale-freeness and heterogeneity of the network by (Aldana
2003). We show that centralized network are capable to produce a number of
patterns, while being protected against environment fluctuations.

Network models usually involve interactions between transcription factors
(TFs) (Reinitz et al. 1991). In the last years, a great attention has been
focused on microRNAs (He and Hannon 2004, Bartel 2009, Hendrikscon et al.
2009, Ihui et al. 2010). MicroRNAs (miRNAs) are short ribonucleic acid (RNA)
molecules, on average only 22 nucleotides long and are found in all eukaryotic
cells. miRNAs are post-transcriptional regulators that bind to complementary
sequences on target messenger RNA transcripts (mRNAs) and repress translation
or trigger mRNA cleavage and degradation. Thus, miRNAs have an impact on gene
expression and it was shown recently that they contribute to canalization of
development (Li et al. 2009). (Shalgi et al. 2007) shows the existence of many
genes submitted to extensive miRNA regulation with many TF among these ”target
hubs”. Without excluding other applications, we consider regulation of TF by
miRNAs as a possible example of a centralized network. For this particular
situation we generalize the TF network models (Reinitz et al. 1991) to take
into account miRNA satellites and the centralized architecture. The
interaction between network nodes is defined by sigmoidal functions that can
be defined by two parameters: the maximum rates of production $r_{i}$, and
sharpness constants $K_{i}$. Other important parameters, play a key role,
namely, the degradation constants $\lambda_{i}>0$ of centers and satellites.

We obtain a fundamental relation between the main network parameters. This
relation ensures maximal robustness of the network with respect to random
internal and external fluctuations, given a certain amount of flexibility
defined as the number of attractors that are accessible to the network
dynamics. Our mathematical results have a transparent biological
interpretation: centralized motifs can be simultaneously flexible and robust.
One can expect that miRNA molecules, being smaller with respect to TF, are
more mobile and react faster to perturbations. This property plays a key role
in the flexible and robust functioning of centralized motifs.

The paper is organized as follows. Centralized networks are introduced in
Section 2. We also formulate here an important assertion on the flexibility of
general centralized networks. We show that these networks are capable to
generate practically all dynamics, chaotic or periodic, with any number of
equilibrium states. To study robustness with respect to random fluctuations,
in Section 3 we consider a toy model of simple centralized TF - miRNA networks
with a single center. We show here, by an elementary way, that centralized
networks with mutually repressive hub-satellite interaction can produce many
different robust patterns.

##  2 Centralized networks

Centralized networks have been empirically identified in molecular biology,
where the centers can be, for example, transcription factors, while the
satellite regulators can be small regulatory molecules such as microRNAs (Li
et al. 2010). Notice that, in the last decades, the theory of so-called scale-
free networks has become very popular. Scale-free networks (Barabasi and
Albert 2002, Lesne 2006) occur in many areas, in economics, biology and
sociology. In the scale-free networks the probability $P(k)$ that a node is
connected with $k$ neighbors, has the asymptotics $Ck^{-\gamma}$, with
$\gamma\in(2,3)$. Such networks typically contain a few strongly connected
nodes and a number of satellite nodes. Hence, scale-free networks are, in a
sense, centralized.

In order to model dynamics of centralized networks we adapt a gene circuit
model proposed to describe early stages of Drosophila (fruit-fly)
morphogenesis (Mjolness et al. 1991, Reinitz and Sharp 1995). To take into
account the two types of the nodes, we use distinct variables $v_{j}$, $u_{i}$
for the centers and the satellites. The real matrix entry $A_{ij}$ defines the
intensity of the action of a center node $j$ on a satellite node $i$. This
action can be either a repression $A_{ij}<0$ or an activation $A_{ij}>0$.
Similarly, the matrices ${\bf B}$ and ${\bf C}$ define the action of the
centers on the satellites and the satellites on the centers, respectively. Let
us assume that a satellite can not act directly on another satellite (the
principle of divide et impera). We also assume that satellites respond more
rapidly to perturbations and are more diffusive/mobile than the centers. Both
these assumptions are natural if we identify satellites as microRNAs.

Let $M,N$ be positive integers, and let ${\bf A},{\bf B}$ and ${\bf C}$ be
matrices of the sizes $N\times M,M\times M$ and $M\times N$ respectively. We
denote by ${\bf A}_{i},{\bf B}_{j}$ and ${\bf C}_{j}$ the rows of these
matrices. To simplify formulas, we use the notation

$\sum_{j=1}^{M}A_{ij}v_{j}={\bf A}_{i}v,\quad\sum_{l=1}^{M}B_{jl}v_{l}={\bf
B}_{j}v,\quad\sum_{k=1}^{N}C_{jk}u_{k}={\bf C}_{j}u.$

Then, the network model reads (we exclude diffusion effects):

$\frac{du_{i}}{dt}=\tilde{r}_{i}\sigma\left({\bf
A}_{i}v-\tilde{h}_{i}\right)-\tilde{\lambda}_{i}u_{i},$ (1)
$\frac{dv_{j}}{dt}=r_{j}\sigma\left({\bf B}_{j}v+{\bf
C}_{j}u-h_{j}\right)-\lambda_{j}v_{j}.$ (2)

We assume that the rate coefficients $r_{j},\tilde{r}_{i}$ are non-negative:
$r_{i},\tilde{r}_{i}\geq 0$. Here $i=1,...,N$, $j=1,...,M$ and $\sigma$ is a
monotone and smooth (at least twice differentiable) sigmoidal function such
that

$\sigma(-\infty)=0,\quad\sigma(+\infty)=1.$ (3)

Typical examples can be given by the Fermi and Hill functions:

$\sigma(x)=\frac{1}{1+\exp(-x)},\quad\sigma_{H}(x)=\frac{x^{p}}{K_{a}^{p}+x^{p}},$
(4)

where $K_{a}$, $p>0$ are parameters and in the second case $x>0$. For $x<0$ we
set $\sigma_{H}(x)=0$. Analytical and computer simulation results are similar
for both variants $\sigma$ and $\sigma_{H}$.

The parameters $\lambda_{i},\tilde{\lambda}_{i}$ are degradation coefficients,
and $h_{i},\tilde{h}_{i}$ are thresholds for activation.

Let us prove that the gene network dynamics defines a dissipative dynamics. In
fact, there exists an absorbing set ${\cal B}$ defined by

${\cal B}=\\{w=(u,v):0\leq v_{j}\leq r_{j}\lambda_{j}^{-1},\ 0\leq
u_{i}\leq\tilde{r}_{i}\tilde{\lambda}_{i}^{-1},\ j=1,...,M,\ i=1,...,N\\}.$

One can show, by comparison principles for ordinary differential equations,
that

$\begin{split}0\leq
u_{i}(x,t)\leq\tilde{\phi}_{i}(x)\exp(-\tilde{\lambda}_{i}t)+\tilde{r}_{i}\tilde{\lambda}_{i}^{-1}(1-\exp(-\tilde{\lambda}_{i}t)),\\\
0\leq
v_{i}(x,t)\leq\phi_{i}(x)\exp(-\lambda_{i}t)+r_{i}\lambda_{i}^{-1}(1-\exp(-\lambda_{i}t)).\end{split}$
(5)

Therefore, solutions of (1), (2) exist for all times $t$ and they enters for
the set ${\cal B}$ at a time moment $t_{0}$ and then stays in this set for all
$t>t_{0}$. So, our system defines a dissipative dynamics and all
concentrations are positive if they are positive at the initial moment. In
mathematical terms, the Cauchy problem (initial value problem) for our system
is well posed.

##  3 Complex dynamics of centralized networks

Let us show that the centralized networks have a formidable power in dynamics
generation. First, we will find an asymptotic simplification of the dynamics,
then show that any dynamics, periodic, chaotic, or with a number of stable
steady states can be approximated by centralized networks.

###  3.1 Simplified dynamics when satellites are fast

We suppose here that the $u$-variables are fast and the $v$-ones are slow.
Then the fast $u$ variables are slaved, for large times, by the slow $v$
modes: one has $u=U(v)+\tilde{u}$, where $\tilde{u}$ is a small correction.
This means that, for large times, the satellite dynamics is defined almost
completely by the center dynamics.

To realize this approach, let us assume that the parameters of the system
satisfy the following conditions:

$\ |A_{jl}|,|B_{il}|,|C_{ij}|,|\tilde{h}_{i}|,|h_{j}|<C_{0},$ (6)

where $i=1,2,...,N,\ \ i,l=1,...,M,\ j=1,...,N$,

$0<C_{1}<\tilde{\lambda}_{j},$ (7)

and

$r_{i}=\kappa R_{i},\quad\tilde{r}_{i}=\kappa\tilde{R}_{i},$ (8)

where

$|R_{i}|,|\tilde{R}_{i}|<C_{5},\quad\lambda_{i}=\kappa\bar{\lambda}_{i},\
|\bar{\lambda}|<C_{6},$ (9)

where $\kappa$ is a small parameter, and where all positive constants $C_{k}$
are independent of $\kappa$.

Assertion 2.1. Under assumptions (6), (7), (8) for sufficiently small
$\kappa<\kappa_{0}$ solutions $(u,v)$ of (1), (2) satisfy

$u=U(v(t))+\tilde{u}(t),$ (10)

where the $j$-th component $U_{j}$ of $U$ is defined by

$-\tilde{\lambda}_{j}U_{j}=\kappa G_{j}(v),$ (11)

where

$G_{j}=\tilde{R}_{j}\sigma\left({\bf A}_{j}v(t)-\tilde{h}_{j}\right)$

The function $\tilde{u}$ satisfies estimates

$|\tilde{u}|<c\kappa^{2}+R\exp(-\beta t),\quad\beta>0.$ (12)

The $v$ dynamics for large times $t>C_{1}|\log\kappa|$ takes the form

$\frac{dv_{i}}{dt}=\kappa F_{i}(u,v)+w_{i},$ (13)

where $w_{i}$ satisfy

$|w_{i}|<c\kappa^{2}$

and

$F_{i}(u,v)=R_{i}\sigma\left({\bf B}_{i}v+{\bf
C}_{i}U(v)-h_{i}\right)-\bar{\lambda}_{i}v_{i}.$

This assertion, known in computational biology as the quasi-steady state
assumption, can be proved by well known methods from the theory of
differential equations (Henry 1981).

###  3.2 Realization of prescribed dynamics by networks

Our next goal is to show that dynamics (13) can realize, in a sense, arbitrary
dynamics of the centers. To precise this, let us describe the method of
realization of the vector fields for dissipative systems (proposed by Poláčik
1991, for applications see, for example, Dancer - Poláčik 1999, Rybakowski
1994, Vakulenko 2000). This method is based on the well developed theory of
invariant and inertial manifolds, see Marion 1989, Mane 1977, Constantin,
Foias, Nicolaenko and Temam, 1989, Chow-Lu 1988, Babin-Vishik 1988). One can
show that there are systems enjoying the following properties:

A These systems generate global semiflows $S_{\cal P}^{t}$ in an ambient phase
space $H$. These semiflows depend on some parameters $\cal P$ (which could be
elements of another parameter space $\cal B$). They have global attractors and
finite dimensional local attracting invariant $C^{1}$ (continuously
differentiable) - manifolds $\cal M$, at least for some $\cal P$.

B Dynamics of $S^{t}_{\cal P}$ reduced on these invariant manifolds is, in a
sense, ”almost completely controllable”. It can be described as follows.
Assume the differential equations

$\frac{dp}{dt}=F(p),\quad F\in C^{1}(B^{n})$ (14)

define a dynamical system in the unit ball ${B}^{n}\subset{\bf R}^{n}$.

For any prescribed dynamics (14) and any $\delta>0$, we can choose suitable
parameters ${\cal P}={\cal P}(n,F,\delta)$ such that

B1 The semiflow $S_{\cal P}^{t}$ has a $C^{1}$\- smooth locally attracting
invariant manifold ${\cal M}_{\cal P}$ diffeomorphic to the ball ${B}^{n}$;

B2 The reduced dynamics $S_{\cal P}^{t}|_{{\cal M}_{\cal P}}$ is defined by
equations

$\frac{dp}{dt}=\tilde{F}(p,{\cal P}),\quad\tilde{F}\in C^{1}(B^{n})$ (15)

where the estimate

$|F-\tilde{F}|_{C^{1}({B}^{n})}<\delta$ (16)

holds. In other words, one can say that, by $\cal P$, the dynamics can be
specified to within an arbitrarily small error.

Thus, all dynamics can occur as inertial forms of these systems. Such systems
can be named maximally dynamically flexible, or, for brevity, MDF systems.

Such dynamics can be chaotic. There is a rather wide broad in different
definitions of ”chaos”. In principle, one can use here any concept of chaos.
If this chaos is stable under small $C^{1}$ -perturbations this kind of chaos
occurs in the dynamics of MDF systems. To fix ideas, we use here, following
Ruelle and Takens 1971, Newhouse, Ruelle and Takens 1971 Smale 1980, Anosov
1995), such a definition. We say that a finite dimensional dynamics is chaotic
if this generates a non-quasiperiodic hyperbolic invariant set $\Gamma$. If,
moreover, this set $\Gamma$ is attracting we say that $\Gamma$ is a chaotic
(strange) attractor. (For definition of hyperbolic sets, see Ruelle 1989,
Anosov 1995). In this paper, we use only the following basic property of
hyperbolic sets, so-called Persistence (Ruelle 1989, Anosov 1995). This means
that the hyperbolic sets are, in a sense, stable(robust): if (14) generates
the hyperbolic set $\Gamma$ and $\delta$ is sufficiently small, then dynamics
(14) also generates another hyperbolic set $\tilde{\Gamma}$. Dynamics (14) and
(15) restricted to $\Gamma$ and $\tilde{\Gamma}$ respectively, are
topologically orbitally equivalent (on definition of this equivalence, see
Ruelle 1989, Anosov 1995). It is important to mention that a chaos in
dissipative systems may be stable, in the sense of structural stability, and
although not yet observed in gene networks, structurally stable chaotic
itineracy is thought to play a functional role in neuroscience (Rabinovitch
1998).

Therefore, any possible chaotic robust dynamics can be generated by the MDF
systems, for example, the Smale horseshoes, Anosov flows, the Ruelle-Takens-
Newhouse chaos, see Newhouse, Ruelle, and Takens, 1971, Smale 1980, Ruelle
1989. Some examples of the MDF systems were given in Dancer- Poláčik 1999,
Rybakowski 1994, Vakulenko 2000.

Assertion 2.1 allows us to apply this approach to centralized network
dynamics. To this end, assume that (8) and (9) hold. Moreover, let us assume

$\lambda_{i}=\kappa^{2}\bar{\lambda}_{i},\quad h_{i}=\kappa\bar{h}_{i}$ (17)

where all coefficients $\bar{h}_{i}$ are uniform in $\kappa$ as $\kappa\to 0$.
We also assume that all direct interactions between centers are absent, ${\bf
B}={\bf 0}$. This constraint is not essential.

Since $U_{j}=O(\kappa)$ for small $\kappa$, we can use the Taylor expansion
for $\sigma$ in (13). Then these equations reduce to

$\frac{dv_{i}(\tau)}{d\tau}=\rho_{i}({\bf
C}_{i}V(v)-\bar{h}_{i})-\bar{\lambda}_{i}v_{i}+\tilde{w}_{i}(t),$ (18)

where $\rho_{i}=\bar{r}_{i}\sigma^{\prime}(0)$, $i=1,2,...,M$ and $\tau$ is a
slow rescaling time: $\tau=\kappa^{2}t$. Due to conditions (17), the
corrections $\tilde{w}_{i}$ satisfy

$|\tilde{w}_{i}|<c\kappa.$

Let us focus now our attention to non-perturbed equation (18) with
$\tilde{w}_{i}=0$. Let us fix the number of centers $M$. The number of
satellites $N$ will be considered as a parameter.

The next important assertion immediately follows from well known approximation
theorems of the multilayered network theory, see, for example, Barron 1993,
Funahashi and Nakamura 1993.

Assertion 2.2. Given a number $\delta>0$, an integer $M$ and a vector field
$F=(F_{1},...,F_{M})$ defined on the ball $B^{M}=\\{|v|\leq 1\\}$, $F_{i}\in
C^{1}(B^{M})$, there are a number $N$, an $N\times M$ matrix ${\bf A}$, an
$M\times N$ matrix ${\bf C}$ and coefficients $h_{i}$, where $i=1,2,...,N$,
such that

$|F_{j}(\cdot)-{\bf C}_{j}W(\cdot)|_{C^{1}(B^{M})}<\delta,$ (19)

where

$W_{i}(v)=\sigma\left({\bf A}_{i}v-h_{i}\right),$ (20)

where $v=(v_{1},...,v_{M})\in{\bf R}^{M}$.

This assertion gives us a tool to control network dynamics. Assume
$\bar{h}_{i}=0$. Then equations (18) with $\tilde{w}_{i}=0$ reduce to the
Hopfield-like equations for variables $v_{i}\equiv v_{i}(\tau)$ that depend
only on $\tau$:

$\frac{dv_{l}}{d\tau}={\bf K}_{l}W(v)-\bar{\lambda}_{l}v_{l},$ (21)

where $l=1,...,M$, the matrix $\bf K$ is defined by
$K_{lj}=\rho_{l}C_{lj}R_{j}\tilde{\lambda}_{j}^{-1}$. The parameters $\cal P$
of (21) are $\bf K$, $M$, $h_{j}$ and $\bar{\lambda}_{j}$.

In this case one can formulate the following result.

Assertion 2.3. Let us consider a $C^{1}$-smooth vector field $Q(p)$ defined on
a ball $B_{R}\subset{\bf R}^{M}$ and directed strictly inside this ball at the
boundary $\partial B^{M}$:

$F(p)\cdot p<0,\quad p\in\partial B^{M}.$ (22)

Then, for each $\delta>0$, there is a choice of parameters $\cal P$ such that
(21) $\delta$ -realizes system (14). This means that (21) is a MDF system.

This follows from Assertions 2.1 and 2.2.

##  4 A toy model of centralized network

In this section we consider a simple centralized network that, nonetheless,
can produce a number of point attractors (stable steady states). Due to its
simple structure, we can investigate here the robustness of this system.

Let us consider a central node interacting with many satellites. This motif
can appear as a subnetwork in a larger scale-free network. In order to study
robustness, we add noise to the model. We consider two types of stochastic
perturbations. The first type of perturbations is a Langevin type additive
noise that can simulate intrinsic stochastic fluctuations of gene expression
dynamics. The choice of additive noise is for the sake of simplicity, however
more general multiplicative noise can be used with no change of the results.
The second type of noise is a shot-like perturbation that can simulate the
external contributions to noise, caused by the environment. Furthermore, we
replace the sigmoid in (2) by a linear function. This is justified in TF -
miRNAs networks, where the action of satellites (miRNA’s) on centers (TF’s) is
post-transcriptional and produces a modulation of the production rate of the
center protein. This modulation can be modeled by a soft sigmoid or even by a
linear function. Moreover, to simplify our model, we assume that all
satellites are, in a sense, equivalent.

The network dynamics can be described then by the following equations:

$\frac{du_{i}}{dt}=-\lambda u_{i}+f_{i}(v)+\xi_{i}(t),$ (23)
$\frac{dv}{dt}=-\nu v+Q(u)+\xi_{0}(t),$ (24)

where $f_{i}$, $Q$ are defined by

$Q(u)=a_{0}+a\sum_{i=1}^{n}u_{i},\quad f_{i}(u)=r\sigma(b(v-h_{i})),$

Here $\xi_{i}$ are noises, the coefficient $\lambda>0$ is a satellite mobility
(degradation rate), $r>0$ is the satellite maximum production rate, $b$
defines a sharpness of center action on the satellites, $\nu>0$ is a center
mobility (degradation rate), $a$ is the strength of the satellites feedback
action on the center.

We consider the following type of noises: non-correlated white noise

$\langle\xi_{i}(t),\xi_{j}(t^{\prime})\rangle=\beta_{i}\delta_{ij}\delta(t-t^{\prime})$
(25)

where $\beta_{i}>0$ are intensities, and shot-like noise

$\xi_{i}(t)=\beta_{i}\eta_{i}\delta(t-\tau_{j})$ (26)

where $\tau_{j}$ are random shot times following a Poisson process,
$\beta_{i}$ are noise amplitude coefficients, and $\eta_{i}$ are random
variables distributed uniformly on $[0,1]$. In numerical simulations we set
$\delta(t-\tau_{j})=1$ with a probability $p_{0}<<1$ and
$\delta(t-\tau_{j})=0$ with the probability $1-p_{0}$, where $\tau_{j}=j\delta
t$, $\delta t$ is a time step. Such noises $\xi_{i}$ can summarize the effect
of a strong environment fluctuations on the satellite and center expression.

We study the problem under the following

Assumption. Let the derivatives of $f_{i}$ and $Q$ satisfy

$f_{i}^{{}^{\prime}}(v)Q^{\prime}(u)>0\ for\ all\ i,u,v.$

Then one can show, following (Hirsch, 1988) that the dynamics is monotone,
and, therefore, all trajectories converge to equilibria. The numerical
simulations confirm this fact. Notice that the above assumption is not needed
when satellites are fast, because in this case the asymptotic dynamics is one
dimensional and in dimension one all the attractors are stable steady states
(point attractors). Although this simple system can not generate chaos or
periodic behavior, the number of point attractors can be arbitrarily large,
and thus this system is nonetheless flexible.

###  4.1 Multistationarity of the toy model

Let us fix the signs of the satellite actions on the center assuming that
$a<0$. This restriction is fulfilled in gene networks, where the centers are
transcription factors (TF) and the satellites are microRNAs (indeed, usually
microRNA can only repress transcription factors). Let us show that the toy
model admits coexistence of any number of point attractors.

Let us make a transformation reducing (23) and (24) to a system of two
equations introducing a new variable $Z$ by

$Z=\sum_{j=1}^{N}u_{j},\quad G(v)=\sum_{j=1}^{N}r\sigma(b(v-h_{j})).$

Then, by summarizing eqs. (23), one obtains

$\frac{dZ}{dt}=-\lambda Z+G(v),$ (27) $\frac{dv}{dt}=-\nu v+a_{0}+aZ.$ (28)

This system is relatively simple and it can be studied analytically and
numerically. Since all trajectories are convergent we obtain that the
attractor consists of equilibria defined by

$\lambda Z+F(Z)=0,\quad F(Z)=G(\nu^{-1}(a_{0}+aZ)).$ (29)

Let ${\cal P}=\\{b,h_{i},N,a_{0},a\\}$ be free parameters that can be
adjusted. Like to the previous section, we can ”control” the nonlinearity $F$
by ${\cal P}$ and use the fact that $F(Z)$ can approximate arbitrary smooth
functions. The following assertion shows that the system is multi-stationarity
with an arbitrary number of point attractors:

Assertion 3.1. Let $N$ be a positive integer. Then there are coefficients
$b,\lambda>0,r>0$, where $i=1,...,N$, $\nu>0$ and $h_{i},a_{0},a$ in such a
way that equation (29) has at least $n+1$ stable roots that can be placed in
any given positions in the $Z$-space.

The main idea of the proof can be illustrated by Fig. 1 and holds in both
cases of the Fermi and the Hill sigmoids. Let us make a variable change
$w=\lambda Z/r$. The steady states are solutions of the equation $rw=F(w)$,
where the function $F(w)$ is close to a step function with $N$ steps; each
step is given by the function $\sigma(\gamma(w-\bar{h}_{i}))$ that is close to
Heaviside step function for large $\gamma$. Here $\gamma$ is a parameter that
defines the sigmoid sharpness:

$\gamma=abr(\nu\lambda)^{-1}.$ (30)

The steady states of the system are given by the intersections between the
graph of $F(w)$ and the straight line of slope $r$. An elementary argument
shows that the intersections lying on horizontal segments of the graph of
$F(w)$ are stable attractors, whereas the intersections on ascending vertical
segments correspond to repellers.

The position of the $i$-th step in $w$-space is $\bar{h}_{i}$ and its height
is $r$. Under an appropriate choice of $\bar{h}_{i}$ this entails our
assertion (see Fig. 1). In the neural network theory, $\gamma$ is known as
gain parameter. This quantity, defined as the product of rates on sharpness
divided on the product of degradation coefficients, gives the maximal possible
density of the equilibrium states in $w$-space.

It is useful to note that one gets $n+1$ attractors on the horizontal segments
of the step function provided that $\bar{h}_{i}$ decrease with $i$.

Notice that the main condition to obtain flexibility (multistationarity) is
the sharpness of the sigmoidal function, meaning that the gain parameter
$\gamma$ should be large. The construction is robust: we can vary
$w_{i},b,\bar{h}_{i}$ but the number of equilibria is conserved.

###  4.2 Robustness and stability of attractors

The roots of Eq.(29) are point attractors and then they are dynamically
stable, otherwise, they are repellers and unstable. In Fig. 1, attractors
correspond to intersections of the straight line $y=rw$ with the curve
$y=F(w)$, lying on horizontal segments of the graph of $F$. A simple argument
suggests that the positions of these attractors are robust with respect to
variations of the thresholds $\bar{h}_{i}$. Indeed, a perturbation of
$\bar{h}_{i}$ induces a horizontal shift of the step
$\sigma(\gamma(w-\bar{h}_{i}))$, and the positions of the attractors are only
slightly affected.

Figure 1: Intersections of the curve $y=F(w)$ and the straight line $y=rw$
correspond to steady states; intersections on horizontal segments of the graph
of $F$ correspond to stable steady states.

More insight into robustness of the centralized toy model can be obtained by
considering the noisy case $\xi_{i}\neq 0$.

First, let us consider the case of the Langevin noise (25). We are interested
in the robustness of the number and positions of the attractors with respect
to noises $\xi_{i}(t)$. Near a point attractor, the equations (27), (28) can
be linearized. The linearized dynamics is defined by the following matrix
${\bf H}$:

$\left(\begin{array}[]{ccc}-\lambda&\mu\\\ a&-\nu\\\ \end{array}\right)$

where $\mu=G^{\prime}(v_{eq})$. For large $\gamma$ and for stable stationary
states $\mu$ is small, $\mu=G^{\prime}(v_{eq})\to 0$ as $\gamma\to\infty$. Let
us assume that the noises $\xi_{i}$ are independent white noises. Using
standard results from the theory of linear stochastic differential equations,
see, for example, Keizer 1987, it follows that small deviations $\delta
Z,\delta v$ from the equilibrium are normally distributed with the density

$\rho(\delta Z,\delta v)=const\exp(-X\cdot{\bf M}^{-1}\cdot X^{tr}),\quad
X=(\delta Z,\delta v),$ (31)

where ${\bf M}$ is a symmetric, positively defined, $2\times 2$ covariation
matrix with entries $m_{11},m_{22},m_{12}=m_{21}$. This matrix can be defined
by the well known relation (the fluctuation-dissipation theorem):

${\bf HM+MH}^{tr}=-{\bf B},$ (32)

where ${\bf B}=diag(B_{1}^{2},B_{2}^{2})$ and, since the noises are non-
correlated, $B_{1}=\sqrt{\sum_{i=1}^{N}\beta_{i}^{2}}$, $B_{2}=\beta_{0}.$

As a result, a characteristic fluctuation amplitude $F_{A}$ is proportional to
the maximum $\max\\{\theta_{1}^{1/2},\theta_{2}^{1/2}\\}$ where $\theta_{i}$
are eigenvalues of ${\bf M}$. Eq. (32) can be resolved explicitly and
$\theta_{i}$ can be found.

Now we can investigate the following problem: how to tune the parameters
$\lambda,\nu$ and $a$ to obtain the minimal fluctuation amplitude $F_{A}$ with
respect to the noise under a given multistationarity level (this means
$\gamma=\gamma_{0}>>1$ is fixed but we can vary the degradation rates
$\nu,\lambda$). This optimization problem can be resolved numerically. The
results, which describe the optimal $\lambda_{opt},\nu_{opt}$ as functions of
$\gamma$, are as follows.

The case A), $B_{2}>>B_{1}$, the center is under a stronger noise than the
satellites. Then the center degradation rate $\nu_{opt}$ should be large, and
$\lambda_{opt}$ is a small, decreasing in $\gamma$ function.

The case B), $B_{2}\leq B_{1}$, the center is under smaller noise than the
satellites. Then the center degradation rate $\nu_{opt}$ should be smaller,
$\lambda_{opt}>\nu_{opt}$, and the both parameters are decreasing in $\gamma$.
This situation is illustrated by Fig. 2.

Figure 2: Optimal degradation parameters $\nu_{opt}$, $\lambda_{opt}$
(minimizing the eigenvalues of the matrix ${\bf H}$) as functions of the gain
parameter $\gamma$ in the case $B_{2}\leq B_{1}$, when the center is under
smaller noise than the satellites.

The classical ideas of the invariant manifold theory, discussed in the
preceding section, allow us to systematize these results. The centralized
network can function under two main and quite opposite regimes. The first one
arises when $\lambda>>\nu$. Then the satellite dynamics is slaved by the
center motion. The center dominates and such a regime can be named power of
the center. This regime is stable if $\beta_{0}$ is small, but $\beta_{i}$ are
large (the noises act on satellites mainly, case B). Considering that the
noise intensity is larger for those components that are expressed in larger
copy numbers, the case should be representative for miRNA-TF networks, when
miRNA are in smaller copy numbers than the transcription factors. In this case
the noise perturb satellite states ($u_{i}$) but, since the satellites are
controlled by the center state $v$, satellites return to the normal states and
dynamics is robust, the noise does not damage the attractor. The opposite
regime is when $\lambda<<\nu$. Then, opposite to the previous situation, the
center dynamics is slaved by the satellites motion. Such a regime can be named
satellite democracy. This regime is stable when $\beta_{0}$ is large, but
$\beta_{i}$ are small (the noise acts stronger on the center, case A). Here
the noise can perturb the center state ($v$) but this state can be restored by
satellites. The large time dynamics is robust, again the noise does not damage
the attractor.

Similar results are illustrated in the case of a shot noise in Fig. 3.

So, we obtain an interesting connection between robustness, multistationarity
and network rate: to support robustness and multistationarity in a noisy
situation, we should decrease the degradation constants. Multistationarity of
molecular switches is important in decision making processes in
differentiation, development, and immune response of the organisms. Our
finding means that noise protected switches are necessarily slow.

Figure 3: Numerical simulations of the system’s trajectories under shot noise.
The parameters were as follows: $N=6$, $\sigma(z)=\sigma_{H}(z)$ with $p=4$,
$b=20$, $K_{a}=1$, $h_{i}=i$, $r=1$, $\lambda=5$, $a_{0}=0.3\nu$, $x\in[0,10]$
and $t\in[0,120]$. The parameter $\beta_{i}=0$ for all $i$ beside $i=3$, where
$\beta_{3}=50$. The two functioning regimes correspond to different values of
$\nu,a$, namely, $a=50,\ \nu=5$ (satellites democracy (SD), a)), and $a=5,\
\nu=0.5$ (power of the center (PC), b)). In the both cases the system shows
multistationarity. For the chosen initial data trajectories converge to an
attractor $v\approx 5.3$ as $t>>1$ in the PC regime and also in the SD regime
$v\to 3.3$. This means that the fast center loses the attractor control, while
the slow center controls dynamics even under large deviations.

##  5 Conclusion

We have considered networks with two types of nodes. The $v$-nodes, called
centers, are hyperconnected and interact one to another via many $u$-nodes,
called satellites. We show, by recently advanced mathematical methods, that
this centralized network architecture, allows us to control network dynamics
to create complicated dynamical regimes. This network organization creates
feedback loops that are capable to generate practically all kinds of dynamics,
chaotic or periodic, or having a number of equilibrium states. This strong
flexibility could be crucial for adaptive biological functions of these
networks.

Using the simple example of a motif with a single center, we also argued that
centralized networks can perform trade-offs between flexibility and
robustness. To support both flexibility and robustness in a noisy situation,
the network should function in a slow manner, i.e, we propose slow-down as a
way to increase stability.

Which of the nodes should be slowed-down depends on the fluctuations. Basic
ideas from the invariant manifold theory show that if the noises act on the
satellites, then, in order to conserve dynamics and the attractor structure,
the center should be slow and controls the satellites (we called this regime
power of the center). In the opposite case, when the noise acts on the center,
the satellites should be slow in order to control the center and the global
dynamics (we called this regime satellites democracy).

We did not consider here extrinsic noise or parametric variability of the
system, that we plan to study in the future. We also think that the slow-down
effect could be observed in all systems where there is a separation into slow
and fast variables, independently of architecture.

Acknowledgements. The authors are grateful to Maria Samsonova and Vitaly
Gursky for useful discussions. We are thankful to M. S. Gelfand and his
colleagues for interesting discussions in Moscow.

The first author was supported by the Russian Foundation for Basic Research
(Grant Nos. 10-01- 00627 s and 10-01-00814 a) and the CDRF NIH (Grant No.
RR07801). We are grateful to the anonymous referees for important remarks.

## References

  * [1]
  * [2] R. Albert and A. L. Barabási (2002),  Rev. Modern Physics, 74, 47-97
  * [3]
  * [4] M. Aldana, (2003),  Boolean dynamics of networks with scale-free topology, Physica D 185, 45-66
  * [5]
  * [6] A. B. Babin and M. I. Vishik,  Regular attractors of semigroups and evolution equations,  J. Math. Pures Appl. . 62 441 – 491, 1983.
  * [7]
  * [8] A. Barron (1993),  Universal Approximation Bounds for superpositions of a sigmoidal functions,  IEEE Trans. on Inf. theory, 39, 930-945. 1993
  * [9]
  * [10] D. P. Bartel,  MicroRNAs: target recognition and regulatory functions  Cell, 136 (2), 215–33.
  * [11]
  * [12] E. N. Dancer and P. Poláčik,  Realization of vector fields and dynamics of spatially homogeneous parabolic equations.  Memoirs of Amer. Math. Society, 140, no. 668, 1999.
  * [13]
  * [14] Dynamical Systems with Hyperbolic Behaviour, D. V. Anosov (ed). (Dynamical Systems 9).  Encyclopedia of Mathematical Sciences Vol. 66. Translated from Russian., Springer V., Berlin, Heidelberg, New-York (1995).
  * [15]
  * [16] C. Foias, G. Sell and Temam R.,  Inertial Manifolds for nonlinear evolutionary equations,  J. Diff. Equations, 73, 309–353, 1988.
  * [17]
  * [18] K. Funahashi and Y. Nakamura.  Approximation of dynamical systems by continuous time recurrent neural networks.  Neural Networks, 6:801–806, 1993.
  * [19]
  * [20] J. K. Hale,  “ Asymptotic behavior of dissipative systems”,  American Mathematical Society, Providence 1988.
  * [21]
  * [22] L. He, G. J. Hannon,  MicroRNAs: small RNAs with a big role in gene regulation,  Nat. Rev. Genet. 5 (7) 522–31, 2004.
  * [23]
  * [24] D. G. Hendrickson, D. J. Hogan, Heather L. McCullogh, J. W. Myers, D. Herschlag, J. E. Ferrell, P. O. Brown,  Concordant Regulation of Translation and mRNA Abundance for Hundreds of Targets of a Human microRNA  PLOS Biology, 7, 2–19, 2009.
  * [25]
  * [26] D. Henry,  Geometric Theory of Semiliniar Parabolic Equations.  Springer, New York, 1981.
  * [27]
  * [28] M. W. Hirsch,  Stability and convergence in strongly monotone dynamical systems J. Reine. Angew. Math. 383, 1– 58, 1988.
  * [29]
  * [30] J. J. Hopfield,  Neural networks and physical systems with emergent collective computational abilities,  Proc. of Natl. Acad. USA, 79, 2554-2558, 1982.
  * [31]
  * [32] M. Ihui, G. Martello and S. Piccolo,  MicroRNA control of signal transduction  Nature Review, Molecular Cell Biology, 11, 252–263, 2010.
  * [33]
  * [34] S. A. Kauffman, (1969)  Metabolic stability and epigenesis in randomly constructed genetic nets.  Journal of Theoretical Biology, 22, 437-467.
  * [35]
  * [36] J. Keizer,  Statistical Thermodynamics of Nonequilibrium processes.  Springer Verlag, New York etc., 1997
  * [37]
  * [38] A. Lesne  Complex networks: from graph theory to biology.  Letters in Math. Phys., 78, 235-262, 2006.
  * [39]
  * [40] Li Li, Jianzhen Xu, Deyin Yang, Xiarong Tan, Hongfei Wang,  Computational approaches for microRNA studies: a review  Mamm. Genome, 21, 1–12, 2010.
  * [41]
  * [42] H. Meinhardt.  Models of biological pattern formation.  Academic Press, London, 1982.
  * [43]
  * [44] E. Mjolness, D.H. Sharp, and J. Reinitz.  A connectionist model of development.  J. Theor. Biol., 152:429–453, 1991.
  * [45]
  * [46] J.D. Murray.  Mathematical Biology.  Springer, New York, 1993.
  * [47]
  * [48] R. Newhouse, D. Ruelle and F. Takens,  Occurence of strange axiom A attractors from quasi periodic flows,  Comm.Math. Phys. 64 , 35–40, 1971.
  * [49]
  * [50] J. Reinitz and D. H. Sharp, Mechanism of formation of eve stripes,  Mechanisms of Development, 49, 133-158, 1995.
  * [51]
  * [52] D. Ruelle.  Elements of differentiable dynamics and bifurcation theory.  Acad. Press, Boston, 1989.
  * [53]
  * [54] D. Ruelle and F. Takens,  On the nature of turbulence,  Comm. Math. Phys, 20, 167 –192, 1971.
  * [55]
  * [56] K. P. Rybakowski,  Realization of arbitrary vector fields on center manifolds of parabolic Dirichlet BVP’s,  J. Differential Equations 114, 199–221, 1994.
  * [57]
  * [58] S. Smale.  On the differential equations of species in competition.  J.Math.Biol., 3:5–7, 1976.
  * [59]
  * [60] S. Vakulenko,  Dissipative systems generating any structurally stable chaos,  Advances in Diff. Equations, 5, 1139-1178, 2000.
  * [61]
  * [62] Zhao J, Yu H, Luo J. H, Cao Z. W, Li Y. X (2006).  Hierarchical modularity of nested bow-ties in metabolic networks.  BMC Bioinformatics 7, p. 386. doi:10.1186/1471-2105-7-386. PMC 1560398. PMID 16916470.
  * [63]