# Logic-dependent emergence of multistability, hysteresis, and biphasic
dynamics in a “minimal” positive feedback network with an autoloop

Akriti Srivastava<EMAIL_ADDRESS>Department of Mathematics and
Statistics, Indian Institute of Technology Kanpur, Kanpur-208016, India
Department of Mathematics, Dayananda Sagar College of Engineering,
Bangalore-560078, India  Mubasher Rashid Corresponding Author: Mubasher Rashid
<EMAIL_ADDRESS>or<EMAIL_ADDRESS>Department of Mathematics
and Statistics, Indian Institute of Technology Kanpur, Kanpur-208016, India

###### Abstract

Cellular decision-making (CDM) is a dynamic phenomenon crucial for development
and diseases. CDM is often, if not invariably, controlled by regulatory
networks defining interactions between genes and transcription factor
proteins. Traditional studies have focussed on molecular switches such as
positive feedback circuits that upon multimerization exhibit at most
bistability. However, higher-order dynamics such as tristability is also
prominent in many biological processes. It is thus imperative to identify a
”minimal” circuit that can alone explain mono, bi, and tristable dynamics. In
this work, we consider a two-component positive feedback network with an
autoloop and explore these regimes of stability for different degrees of
multimerization and the choice of Boolean logic functions. We report that this
network can exhibit numerous dynamical scenarios such as bi-and tristability,
hysteresis, and biphasic kinetics, explaining the possibilities of abrupt cell
state transitions and the $smooth\ state\ swap$ without a step-like switch.
Specifically, while with monomeric regulation and competitive OR logic, the
circuit exhibits mono-and bistability and biphasic dynamics, with non-
competitive AND and $\mathcal{OR}$ logics only monostability can be achieved.
To obtain bistability in the latter cases, we show that the autoloop must have
(at least) dimeric regulation. In pursuit of higher-order stability, we show
that tristability occurs with higher degrees of multimerization and with non-
competitive $\mathcal{OR}$ logic only. Our results, backed by rigorous
analytical calculations and numerical examples, thus explain the association
between multistability, multimerization, and logic in this ”minimal” circuit.
Since this circuit underlies various biological processes, including
epithelial-mesenchymal transition which often drives carcinoma metastasis,
these results can thus offer crucial inputs to control cell state transition
by manipulating dimerization and the logic of regulation in cells.

Keywords: Cell fate decision; Multistability; Minimal genetic circuit;
Positive feedback loop; Biphasic dynamics; Boolean logic.

##  1 Introduction

Cellular decision-making is a cell non-autonomous process where at each cell-
lineage branch point, a cell drives into one of the alternative distinct cell
types [1]. Investigating dynamical principles of regulatory networks that
govern cellular decision-making is essential to understanding and controlling
stepwise lineage decisions of cells [2, 3, 4]. A key aspect of this decision-
making process is the ability of cells to manifest multiple stable states or
phenotypes in response to varying internal and external cues without altering
their genetic makeup [5, 1]. This phenomenon, known as multi-stability, is
pivotal to cellular differentiation [5, 6], reprogramming [2, 1], and its
manifestation through epithelial-mesenchymal plasticity [7], a cellular
program enabling bidirectional transition among epithelial, mesenchymal, and
one or more hybrid epithelial-mesenchymal cells, enables intra-tumor
heterogeneity [8, 9], induces cancer progression and metastasis of carcinomas
[10, 11]. Since cellular decision-making is largely controlled by regulatory
networks defining “molecular switches” [12, 13], unraveling the dynamical
behaviors and logic of multi-stable switches has profound implications in
synthetic biology and regenerative medicine.

A necessary condition for nonlinear control systems and networks to exhibit
multistability is the presence of positive feedback loops [14, 15, 16, 17,
13]. Interestingly, this requirement also extends to biological systems where
a commonly observed network structure underlying multistability is the toggle
switch [18], comprising mutual inhibition of two opposite fate-determining
transcription factors thus forming a positive feedback loop. This mutual
repression enables cells to adopt different states, driving an ’either-or’
binary choice between alternative cell fates from a common progenitor [2].
However, it is rational to posit that cellular decision-making doesn’t rely
solely on binary outcomes. The intricate structure of gene regulatory networks
enables them to exhibit multiple alternative stable states. Examples of, for
instance, bi-and tristable solutions are found across biological contexts
including development [19], differentiation [20, 21, 22, 23], epigenetic
processes [24], and metastasis of carcinomas [25, 26, 27]. Often, a two-
component system operates at the core of the underlying bistable decision-
making circuits [2, 28, 29]. It is imperative to identify the ”minimal” (in
terms of the number of variables and the regulatory interactions) two-
component genetic circuits capable of demonstrating bistability. Two such
networks are (i) circuits featuring multimeric regulation and lacking
autoloop, and (ii) circuits employing monomeric regulators along with a
singular monomeric autoloop [30, 12]. Jules et al., [31] and Zhang et al.,
[32] analytically proved that such systems can exhibit at most bistability.
This begs a question: are low-dimensional networks such as two-component
circuits capable of only bistable behavior or higher-order stability such as
tristability (three states) can also be possible in these circuits? It may be
noted that compared to bistable dynamics, networks exhibiting tristable
dynamics are not much studied due to the complexity of analytical calculations
as well as numerical simulations.

Recent theoretical and experimental studies found that miR200-ZEB positive
feedback circuit with one multimeric autoloop, which underlies epithelial-
mesenchymal transition, exhibits tristability [27, 26, 25, 33]. It allows
reversible transition among epithelial, mesenchymal, and hybrid epithelial-
mesenchymal cells and is a crucial driver of invasion and metastasis of
carcinomas [34]. However, these studies are purely based on a numerical
approach and lack mathematical rigor. Also, they don’t explore any association
of multimerization with tristability. Motivated by this, we asked a general
question: what can be the minimal ”multimeric” requirements for a two-
component feedback network with one autoloop (Fig. 1a) to exhibit mono-and
bistability and under what conditions can it exhibit tristability? Further,
node $a$ is co-regulated by node $b$ and node $a$ itself through autoloop.
When a gene is regulated by more than one transcription factor, the
combinatorial rule (Boolean logic) must also be considered. We considered
three types of regulatory logic: (a) Competitive OR logic, (b) non-competitive
AND logic, and (c) non-competitive OR logic, and explored their effect on the
order of stability. For simplicity, throughout the text, we denote these three
logics as Type-I OR, AND, and Type-II $\mathcal{OR}$, respectively.

In this work, we thus analyze the dynamics of such a generic, minimal network
(depicted in Fig. 1a) by systematically finding the analytical conditions for
mon, bi, and tristable dynamical regimes. We model the network using a system
of coupled ordinary differential equations for which each stable equilibrium
solution corresponds to a possible stable state (or phenotype or attractor).
We explore how these solutions are affected by the degree of multimerization
as well as the type of regulatory logic. We integrated the autoloop with the
core positive feedback network using three regulatory logics: AND gate, Type-I
OR gate, and Type-II $\mathcal{OR}$ gate, each having a different mathematical
formalism and qualitative behavior (Fig. 1b). Starting with monomeric
regulation, we systematically tested higher degrees of multimerization in each
model. We report that multistability in this network is logic-dependent. In
monomeric regulation with monomeric autoloop, while type-I OR logic exhibits
at most bistability, AND and type-II $\mathcal{OR}$ logic can exhibit only
monostability. To obtain bistability with AND and Type-II $\mathcal{OR}$
logics, we show that a dimeric autoloop is required. Upon testing different
combinations of multimerization, we show that tristability can occur in higher
degrees of multimerization, though bistability is persistently prevalent
irrespective of the logic of regulation. Our results thus comprehensively
provide an association between multistability, multimerization, and logic in
this circuit. Since this circuit underlies various biological processes
including EMT in cancers and the multistability can be mapped to phenotypic
heterogeneity, these results can thus offer crucial inputs to control
phenotypic heterogeneity by manipulating dimerization and synthetically
(re)engineering the logic of regulation in cells.

Figure 1: Illustration of auto-regulated positive feedback network and gate
function response curves. (a) Genes $a$ and $b$ mutually inhibit each other to
form a positive feedback circuit. Also, gene $a$ can promote its expression
(directly or indirectly), which is represented by autoloop. An example of a
biological network corresponding to this schematic diagram is the miR200-ZEB
mutual inhibition circuit underlying EMT and driving metastasis of carcinomas
[27, 26]. (b) Qualitative behavior of three types of Boolean logic functions
modeling co-regulation of gene $a$ (activation by $a$ and inhibition by $b$)
for different values of ($n_{aa}$, $n_{b}$). Type-I OR logic shows mixed
increasing/decreasing responses for different multimerization conditions
(i.e., Hill coefficients; $n_{aa}$, $n_{b}$). The AND and type-II
$\mathcal{OR}$ logic functions show a unified response independent of
multimerization. Also, it is evident that type-I OR and type-II $\mathcal{OR}$
logic functions differ qualitatively. The former in most cases is increasing
and saturating while the latter is always monotonically decreasing.

##  2 Organization of the paper

This paper contains three main sections based on the three different types of
Boolean logic used to model the co-regulation of gene $a$ by two transcription
factors. In section 3, this co-regulation is modeled using competitive OR
logic. In section 4, the non-competitive AND logic is employed. In section 5,
a non-competitive $\mathcal{OR}$ logic is considered to model the co-
regulation. For brevity, throughout the text, competitive OR logic, non-
competitive AND logic, and non-competitive OR logic are referred to as Type-I
OR logic, AND logic, and Type-II $\mathcal{OR}$ logic, respectively. In each
section, we explore the minimal multimerization conditions (manifested through
Hill coefficient values) for the occurrence of mono, bi-and multistability. We
have applied the nullcline approach to show the existence of steady states.
For each model, we analyze the associated characteristic equations and obtain
the conditions for the stability of the steady states and the existence of
saddle-node bifurcations. Numerical simulations are then performed to
illustrate the results. Through bifurcations, we also present possible ways of
state transitions such as step-like abrupt state transition (a manifestation
of saddle-node bifurcation) and biphasic $smooth\ state\ swaps$ (manifested
through sigmoidal curves). In Section 6, we discuss the impact of basal
production rate (leakage) on each model’s dynamics. Finally, section 7
presents a detailed discussion and concluding remarks based on our findings.

##  3 Mathematical Model [Type-I OR Logic]

The model under investigation in this paper is depicted in Fig. 1. The
complete nonlinear mathematical model uses four state variables for the
concentrations of mRNAs and transcription factor proteins (TF). The
concentration of mRNA produced by genes $a$ and $b$ is denoted by $r_{a}$ and
$r_{b},$ respectively, while the corresponding TF protein concentrations are
denoted by $p_{a}$ and $p_{b}.$ The inhibition of gene $a$ and $b$ by TF
$p_{b}$ and $p_{a}$, respectively, is modelled by the Hill function. Also, the
activation of gene $a$ by its own TF, $p_{a}$, is modeled by the Hill function
with a different Hill coefficient, $n_{aa}$. The co-regulation of gene $a$ by
TF $p_{b}$ and also by its own TF $p_{a}$ is modeled using competitive OR
logic [35]. The transcription of genes $a$ and $b$ into mRNA $r_{a}$ and
$r_{b}$ is modeled using the non-linear Hill function while the translation of
these mRNA’s into TF proteins $p_{a}$ and $p_{b}$ is modeled using linear
mass-action kinetics.

Using the above formalism, the following system of coupled nonlinear ordinary
differential equations describe the dynamics of the network in Fig. 1.

$\displaystyle\frac{dr_{a}}{dt}=m_{a}\frac{(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}}{1+(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}+(\frac{p_{b}}{\theta_{b}})^{n_{b}}}-\gamma_{a}r_{a}+A_{1},$
(1)
$\displaystyle\frac{dr_{b}}{dt}=m_{b}\frac{1}{1+(\frac{p_{a}}{\theta_{a}})^{n_{a}}}-\gamma_{b}r_{b}+B_{1},$
$\displaystyle\frac{dp_{a}}{dt}=k_{a}r_{a}-\delta_{a}p_{a},$
$\displaystyle\frac{dp_{b}}{dt}=k_{b}r_{b}-\delta_{b}p_{b},$

with $r_{a}(0)\geq 0,\,r_{b}(0)\geq 0,\,p_{a}(0)\geq 0,\,p_{b}\geq 0,$ where
$m_{i}$ are the maximum transcription rates, $\gamma_{i}$ are the mRNA
degradation rates, $k_{i}$ are the translation rates, and $\delta_{i}$ are the
protein degradation rates for genes $i=a,\,b$. $A_{1}$ and $B_{1}$ represent
the basal production rates (in the absence of activation and inhibition) of
genes $a$ and $b,$ respectively. The threshold parameters
$\theta_{a},\,\theta_{b},\,\theta_{aa}$ denote the thresholds of $p_{a}$ and
$p_{b}$ to induce a significant response of $r_{a}$ and $r_{b}$. The integer
parameters $n_{a},\,n_{b}$ and $n_{aa}$ are Hill coefficients that determine
the steepness of Hill curves.

###  3.1 Mathematical Analysis

This section presents the existence and stability analysis of the steady
states for the model system (1). We use the Routh-Hurwitz criterion [36] to
prove the local asymptotic stability of the steady states, while Sotomayor’s
theorem [36] is used to derive the transversality conditions of saddle-node
bifurcation.

We obtained the steady state $S^{*}=(r_{a}^{*},r_{b}^{*},p_{a}^{*},p_{b}^{*})$
by setting the right hand side of the model (1) to zero. From third and fourth
equations of model (1), we obtained
${r_{a}}^{*}=\frac{\delta_{a}{p_{a}}^{*}}{k_{a}}$ and
${r_{b}}^{*}=\frac{\delta_{b}{p_{b}}^{*}}{k_{b}},$ respectively, as
$\frac{dp_{a}}{dt}=0$ and $\frac{dp_{b}}{dt}=0.$ Using $\frac{dr_{b}}{dt}=0$
and then putting the value of $r_{b},$ we get
${p_{b}}^{*}=\frac{1}{k_{2}}\bigg{\\{}B_{1}+m_{b}\frac{\theta_{a}^{n_{a}}}{\theta_{a}^{n_{a}}+{{p_{a}}^{*}}^{n_{a}}}\bigg{\\}}.$
Finally, using the value of $p_{b}$ and $r_{a}$ in $\frac{dr_{a}}{dt}=0,$ we
get a polynomial in $p_{a}$ as

$\begin{split}\bigg{(}m_{a}p_{a}^{n_{aa}}+(A_{1}-k_{1}p_{a})(\theta_{aa}^{n_{aa}}+p_{a}^{n_{aa}})\bigg{)}\theta_{b}^{n_{b}}k_{b}^{n_{b}}(\theta_{a}^{n_{a}}+p_{a}^{n_{a}})^{n_{b}}+(A_{1}-k_{1}p_{a})&\theta_{aa}^{n_{aa}}\bigg{(}B_{1}(\theta_{a}^{n_{a}}+p_{a}^{n_{a}})\\\
+m_{b}\theta_{a}^{n_{a}}\bigg{)}^{n_{b}}=0,\end{split}$ (2)

where $k_{1}=\frac{\gamma_{a}\delta_{a}}{k_{a}}$ and
$k_{2}=\frac{\gamma_{b}\delta_{b}}{k_{b}}.$ Here, $p_{a}^{*}$ is a positive
real root of (2) which we discuss in different cases explicitly.
Case 1: When $n_{aa}=n_{a}=n_{b}=1,$ then from (1), we get a cubic polynomial
in $p_{a}$

$f(p_{a})=Ap_{a}^{3}+Bp_{a}^{2}+Cp_{a}+D,$ (3)

where
$A=k_{1}\,k_{2}\,\theta_{b},\,B=B_{1}\,k_{1}\,\theta_{aa}+k_{1}\,k_{2}\,\theta_{a}\,\theta_{b}-k_{2}\,\theta_{b}\,\left(A_{1}+\mathrm{m_{a}}-k_{1}\,\theta_{aa}\right)\\\
C=k_{1}\,\theta_{aa}\,\theta_{a}\,\left(B_{1}+\mathrm{m_{b}}\right)-k_{2}\,\left(A_{1}+\mathrm{m_{a}}-k_{1}\,\theta_{aa}\right)\theta_{a}\theta_{b}-A_{1}\,B_{1}\,\theta_{aa}-A_{1}\,k_{2}\,\theta_{b}\,\theta_{aa}$
and
$D=-A_{1}\,\theta_{aa}\,\theta_{a}\,\left(B_{1}+\mathrm{m_{b}}+k_{2}\,\theta_{b}\right).$
The nature of the curve $f$ can be observe from (3) as

  1. (i)

$f(0)<0$ because $D<0,$

  2. (ii)

$f\to\infty\,(-\infty)$ as $t\to\infty\,(-\infty)$ because $A>0,$

  3. (iii)

The number of positive real roots of equation (3) can be seen using Descarte’s
rule of sign, mentioned in below Table (1).

Coefficients |  Number of possible positive real roots of $f(p_{a})=0$
---|---
A | B | C | D
+ | - | - | - | 1
+ | + | - | - | 1
+ | + | + | - | 1
+ | - | + | - | 3,1
Table 1: Number of possible positive real roots of $f(p_{a})=0.$

Therefore, $f(p_{a})=0$ will have either unique or three positive real roots.
Case 2: When $n_{a}=n_{b}=1,\,n_{aa}=2,$ then from (1), we get a quartic
polynomial in $p_{a}$

$f(p_{a})=Ap_{a}^{4}+Bp_{a}^{3}+Cp_{a}^{2}+Dp_{a}+E,$ (4)

where
$A=k_{1}\,k_{2}^{n_{b}}\,\theta_{b}^{n_{b}},\,B=-(m_{a}+A_{1}-k_{1}\theta_{a}^{n_{a}})k_{2}^{n_{b}}\theta_{b}^{n_{b}},\,C=k_{1}B_{1}\theta_{aa}^{n_{aa}}-k_{2}^{n_{b}}\theta_{b}^{n_{b}}(\theta_{a}^{n_{a}}(m_{a}+A_{1})-\theta_{aa}^{n_{aa}}k_{1}),\,D=-k_{2}^{n_{b}}\theta_{b}^{n_{b}}\theta_{aa}^{n_{aa}}(A_{1}-k_{1}\theta_{a}^{n_{a}})-\theta_{aa}^{n_{aa}}(A_{1}B_{1}-k_{1}(B_{1}+m_{b})\theta_{a}^{n_{a}}),$
and
$E=-A_{1}\theta_{aa}^{n_{aa}}\theta_{a}^{n_{a}}(k_{2}^{n_{b}}\theta_{b}^{n_{b}}+B_{1}+m_{b}).$
The number of positive real roots of equation (4) can be seen using Descarte’s
rule of sign, mentioned in below Table (2).

Coefficients |  Number of possible positive real roots of $f(p_{a})=0$
---|---
A | B | C | D | E
+ | - | + | \- or + | - | 3, 1
+ | \+ or - | - | + | - | 3, 1
+ | \+ or - | - | - | - | 1
+ | + | + | \- or + | - | 1
Table 2: Number of possible positive real roots of $f(p_{a})=0.$

Case 3: When $n_{a}=2,\,n_{b}=1,\,n_{aa}=2$ then from (1), we get a quintic
polynomial in $p_{a}$

$f(p_{a})=Ap_{a}^{5}+Bp_{a}^{4}+Cp_{a}^{3}+Dp_{a}^{2}+Ep_{a}+F,$ (5)

where
$A=k_{1}\,{k_{2}}^{n_{b}}\,{\theta_{b}}^{n_{b}},\,B=-{k_{2}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\,\left(A_{1}+\mathrm{m}_{a}\right),\,C=B_{1}\,k_{1}\,{\theta_{aa}}^{n_{aa}}+k_{1}\,{k_{2}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\,\\\
\left({\theta_{a}}^{n_{a}}+{\theta_{aa}}^{n_{aa}}\right),\,D=-A_{1}\,B_{1}\,{\theta_{aa}}^{n_{aa}}-{k_{2}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\,\left({\theta_{a}}^{n_{a}}\,\left(A_{1}+\mathrm{m}_{a}\right)+A_{1}\,{\theta_{aa}}^{n_{aa}}\right),\\\
E=k_{a}\,{\theta_{a}}^{n_{a}}\,{\theta_{aa}}^{n_{aa}}\,\left(B_{1}+\mathrm{m}_{b}+{k_{2}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\right)$
and
$F=-A_{1}\,{\theta_{a}}^{n_{a}}\,{\theta_{aa}}^{n_{aa}}\,\left(B_{1}+\mathrm{m}_{b}+{k_{2}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\right).$
The nature of curve $f$ can be observe from (3) as

  1. (i)

$f(0)<0$ because $F<0,$

  2. (ii)

$f\to\infty\,(-\infty)$ as $t\to\infty\,(-\infty)$ because $A>0,$

  3. (iii)

$f(p_{a})=0$ will have either unique or three or five positive real roots.

To deal with the fifth-degree polynomial with parametric coefficients is quite
difficult, so we use the nullcline method to further analyze the existence of
possible positive real roots. From the system (1), we have

$\displaystyle F_{1}(p_{a},p_{b}):=\frac{m_{a}\
p_{a}^{2}\theta_{b}}{\theta_{aa}^{2}\theta_{b}+p_{a}^{2}\theta_{b}+p_{b}\theta_{aa}^{2}}+A_{1}-\frac{\gamma_{a}\delta_{a}p_{a}}{k_{a}}=0$
(6) $\displaystyle
F_{2}(p_{a},p_{b}):=\frac{m_{b}\theta_{a}^{2}}{\theta_{a}^{2}+p_{a}^{2}}+B_{1}-\frac{\gamma_{b}\delta_{b}p_{b}}{k_{b}}=0,$

First to observe the nature of curve $F_{2}(p_{a},p_{b})=0,$ we reduce it in
the form of

$p_{b}=g_{1}(p_{a}):=\frac{k_{b}}{\gamma_{b}\delta_{b}}\Big{\\{}\frac{m_{b}\theta_{a}^{2}}{(\theta_{a}^{2}+p_{a}^{2})}+B_{1}\Big{\\}}$
($\mathscr{A}$)

The following conclusion can be made from ($\mathscr{A}$)

  1. $\mathscr{A}$(i)

$g_{1}(0)=\frac{k_{b}}{\gamma_{b}\delta_{b}}(m_{b}+B_{1})>\frac{k_{b}}{\gamma_{b}\delta_{b}}B_{1}>0$

  2. $\mathscr{A}$(ii)

$g_{1}(p_{a})\to\frac{k_{b}}{\gamma_{b}\delta_{b}}B_{1}$ as
$p_{a}\to\pm\infty$

  3. $\mathscr{A}$(iii)

$g^{\prime}(p_{a})_{1}<0$, i.e, $g_{1}(p_{a})$ is a decreasing function of
$p_{a}.$

The graph of $F_{2}(p_{a},p_{b})=0$ (or $p_{b}=g_{1}(p_{a})$) is shown in Fig.
2 (a).
Now, to observe the nature of curve $F_{1}(p_{a},p_{b})=0,$ we reduce it in
the form of

$p_{b}=g_{2}(p_{a}):=\frac{\left(-\gamma_{a}\delta_{a}p_{a}^{3}+\left(m_{a}k_{a}+A_{1}k_{a}\right)p_{a}^{2}-\gamma_{a}\delta_{a}\theta_{aa}^{2}p_{a}+A_{1}k_{a}\theta_{aa}^{2}\right)\theta_{b}}{\theta_{aa}^{2}\left(p_{a}-\frac{A_{1}k_{a}}{\delta_{a}\gamma_{a}}\right)\delta_{a}\gamma_{a}}.$
($\mathscr{B}$)

The following conclusion can be made from ($\mathscr{B}$)

  1. $\mathscr{B}$(i)

$g_{2}(p_{a})$ is a rational function and it has an asymptote at
$(p_{a}-\frac{A_{1}^{\prime}k_{a}}{\delta_{a}\gamma_{a}})$.

  2. $\mathscr{B}$(ii)

The positive $p_{a}-$intercept can be find with $g_{2}(p_{a})=0,$ which
implies that the positive real roots of
$-\gamma_{a}\delta_{a}p_{a}^{3}+\left(m_{a}k_{a}+A_{1}k_{a}\right)p_{a}^{2}-\gamma_{a}\delta_{a}\theta_{aa}^{2}p_{a}+A_{1}k_{a}\theta_{aa}^{2}=0$
will provide the positive $p_{a}-$intercept of $g_{1}(p_{a}).$

  3. $\mathscr{B}$(iii)

The $p_{b}$-intercept of $g_{2}(p_{a})$ is negative.

  4. $\mathscr{B}$(iv)

We note that for $p_{a}\in(0,\frac{A_{1}k_{a}}{\delta_{a}\gamma_{a}})$, we
have $g_{2}(p_{a})<0$ so if positive real roots of $g_{2}(p_{a})=0$ exits then
it always lie on the right-hand side of asymptote.
The possibility of unique positive $p_{a}$-intercept of $g_{2}(p_{a})$ is
shown in Fig. 2 (b) and Fig. 2 (c). With unique positive $p_{a}$-intercept,
there can be unique and three positive steady states, as shown in Fig. 3 (a)
and Fig. 3 (b), respectively. In Fig. 2 (d), we have shown the three positive
$p_{a}$-intercept of $g_{2}(p_{a})$, and for this case, we obtain three and
unique positive steady states, as shown in Fig. 3 (c) and Fig. 3 (d),
respectively. Therefore, from the above discussion, we discard the case of
five positive steady states.

Figure 2: For $n_{a}=n_{aa}=1=2,\,n_{b}=1;$ (a) Plot of $g_{1}(p_{a})$ as a
decreasing function of $p_{a}$ (b) The graph depicting unique positive
$p_{a}$-intercept of $g_{2}(p_{a})$ (c) The graph depicting another
possibility of unique positive $p_{a}$-intercept of $g_{2}(p_{a})$ (d) The
graph depicting three positive $p_{a}$-intercept of $g_{2}(p_{a}).$

Figure 3: For $n_{a}=n_{aa}=1=2,\,n_{b}=1;$ (a) Existence of unique positive
steady state corresponding to Fig. 2 (b). (b) Existence of three positive
steady states corresponding to Fig. 2 (c) Existence of three positive steady
states corresponding to Fig. 2 (d) Existence of unique positive steady state
corresponding to Fig. 2 (d).

####  3.1.1 Local stability of steady state

######  Theorem 1.

The steady state $S_{1}^{*}$ is locally asymptotically stable if
$\epsilon_{1}>0,\,\epsilon_{3}>0,\,\epsilon_{4}>0,$ and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}>0,$
where $\epsilon_{i},\,i=1,..,4$ are provided in the proof.

###### Proof.

At $S_{1}^{*},$ the Jacobian matrix of the model system (1) is given by

$J_{S_{1}^{*}}=\left[\begin{array}[]{cccc}-\gamma_{a}&0&X&Y\\\
0&-\gamma_{b}&Z&0\\\ k_{a}&0&-\delta_{a}&0\\\
0&-k_{b}&0&-\delta_{b}\end{array}\right].$

where
$X=\frac{\mathrm{m}_{a}\,n_{aa}\,{\mathrm{p_{a}}}^{n_{aa}-1}\,{\theta_{b}}^{n_{b}}\,{\theta_{aa}}^{n_{aa}}\,\left({\mathrm{p_{b}}}^{n_{b}}+{\theta_{b}}^{n_{b}}\right)}{{\left({\mathrm{p_{a}}}^{n_{aa}}\,{\theta_{b}}^{n_{b}}+{\mathrm{p_{b}}}^{n_{b}}\,{\theta_{aa}}^{n_{aa}}+{\theta_{b}}^{n_{b}}\,{\theta_{aa}}^{n_{aa}}\right)}^{2}},\,Y=-\frac{\mathrm{m}_{a}\,n_{b}\,{\mathrm{p_{a}}}^{n_{aa}}\,{\mathrm{p_{b}}}^{n_{b}-1}\,{\theta_{b}}^{n_{b}}\,{\theta_{aa}}^{n_{aa}}}{{\left({\mathrm{p_{a}}}^{n_{aa}}\,{\theta_{b}}^{n_{b}}+{\mathrm{p_{b}}}^{n_{b}}\,{\theta_{aa}}^{n_{aa}}+{\theta_{b}}^{n_{b}}\,{\theta_{aa}}^{n_{aa}}\right)}^{2}},$
and
$Z=-\frac{\mathrm{m}_{b}\,n_{a}\,{\mathrm{p_{a}}}^{n_{a}-1}\,{\theta_{a}}^{n_{a}}}{{\left({\mathrm{p_{a}}}^{n_{a}}+{\theta_{a}}^{n_{a}}\right)}^{2}}.$
The characteristic equation of $J_{S_{1}^{*}}$ is given as

$\lambda^{4}+\epsilon_{1}\lambda^{3}+\epsilon_{2}\lambda^{2}+\epsilon_{3}\lambda+\epsilon_{4}=0.$
(7)

Here,
$\epsilon_{1}=\gamma_{a}+\gamma_{b}+\delta_{a}+\delta_{b},\epsilon_{2}=\gamma_{a}\gamma_{b}+(\gamma_{a}+\gamma_{b})(\delta_{a}+\delta_{b})+\delta_{a}\delta_{b}-k_{a}X,\,\epsilon_{3}=\delta_{a}\delta_{b}(\gamma_{a}+\gamma_{b})+(\delta_{a}+\delta_{b})(\gamma_{a}\gamma_{b})-k_{a}(\delta_{b}+\gamma_{b})X,$
and $\epsilon_{4}=\gamma_{a}\gamma_{b}\delta_{a}\delta_{b}-k_{a}k_{b}YZ-
k_{a}\gamma_{b}\delta_{b}X.$ By the Routh-Hurwitz criteria, all the roots of
(7) have negative real parts provided
$\epsilon_{1}>0,\,\epsilon_{3}>0,\,\epsilon_{4}>0,$ and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}>0.$
∎

####  3.1.2 Saddle Node Bifurcation

We derive the transversality condition for saddle-node bifurcation considering
$\theta_{b}$ as a bifurcation parameter for system (1). Let $\theta_{b}^{SN}$
is the threshold value of $\theta_{b}$ for the occurrence of saddle-node
bifurcation.The system (1) has two steady states ($S_{2}^{*},\,S_{3}^{*}$) for
$\theta_{b}>\theta_{b}^{SN}(\theta_{b}<\theta_{b}^{SN}),$ collide
($S_{2}^{*}=S_{3}^{*}$) at $\theta_{b}=\theta_{b}^{SN}$ and then disappear
when $\theta_{b}<\theta_{b}^{SN}(\theta_{b}>\theta_{b}^{SN})$. The point of
collision is represented as
$S_{SN}^{*}=(r_{a}^{{SN}^{*}},r_{b}^{{SN}^{*}},p_{a}^{{SN}^{*}},p_{b}^{{SN}^{*}}).$
Define $f=(f_{1},f_{2},f_{3},f_{4})^{T},$ where $f_{1},\,f_{2},\,f_{3},$ and
$f_{4}$ are defined as

$\displaystyle
f_{1}=m_{a}\frac{(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}}{1+(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}+(\frac{p_{b}}{\theta_{b}})^{n_{b}}}-\gamma_{a}r_{a}+A_{1},$
$\displaystyle
f_{2}=m_{b}\frac{1}{1+(\frac{p_{a}}{\theta_{a}})^{n_{a}}}-\gamma_{b}r_{b}+B_{1},$
$\displaystyle f_{3}=k_{a}r_{a}-\delta_{a}p_{a},$ $\displaystyle
f_{4}=k_{b}r_{b}-\delta_{b}p_{b},$

Let the Jacobian matrix $J=Df(S_{SN}^{*},\theta_{b}^{SN})$ has a simple
eigenvalue $\lambda=0$ with eigenvector $v=(v_{1},v_{2},v_{3},v_{4})^{T}$, and
the transpose of the Jacobian matrix $J^{T}$ has an eigenvector
$w=(w_{1},w_{2},w_{3},w_{4})^{T}$ to the eigenvalue $\lambda=0.$ Then the
model system (1) experiences a saddle-node bifurcation at
$(S_{SN}^{*},\theta_{b}^{SN})$, if the following transversality conditions
[36] are satisfied:

$w^{T}f_{\theta_{b}}(S_{SN}^{*},\theta_{b}^{SN})=w_{1}\frac{m_{a}n_{b}\theta_{aa}^{n_{aa}}\theta_{b}^{n_{b}-1}{p_{b}^{n_{b}}}^{{SN}^{*}}{p_{a}^{n_{aa}}}^{{SN}^{*}}}{(\mathcal{X})^{2}}\neq
0.$ ($\mathcal{A}_{1}$)

and

$\displaystyle w^{T}[D^{2}f(S_{SN}^{*},\theta_{b}^{SN})(v,v)]=$
$\displaystyle\frac{-w_{1}m_{a}\theta_{aa}^{n_{aa}}\theta_{b}^{n_{b}}}{\mathcal{X}^{3}}\bigg{\\{}2v_{3}^{2}n_{aa}\theta_{b}^{n_{b}}{p_{a}^{n_{aa}-1}}^{{SN}^{*}}(\theta_{b}^{n_{b}}+{p_{b}^{n_{b}}}^{{SN}^{*}})+v_{4}^{2}n_{b}{p_{a}^{n_{aa}}}^{{SN}^{*}}$
$\displaystyle\Big{(}\theta_{b}^{n_{b}}(\theta_{aa}^{n_{aa}}+{p_{a}^{n_{aa}}}^{{SN}^{*}}){p_{b}^{n_{b}-2}}^{{SN}^{*}}(n_{b}-1)-\theta_{aa}^{n_{aa}}{p_{b}^{2n_{b}-2}}^{{SN}^{*}}(n_{b}+1)\Big{)}$
$\displaystyle-2v_{3}v_{4}n_{b}{p_{b}^{n_{b}-1}}^{{SN}^{*}}\Big{(}\theta_{b}^{n_{b}}{p_{a}^{n_{aa}}}^{{SN}^{*}}-\big{(}\theta_{b}^{n_{b}}+{p_{b}^{n_{b}}}^{{SN}^{*}}\big{)}\theta_{aa}^{n_{aa}}\Big{)}\bigg{\\}}$
$\displaystyle-\frac{w_{2}v_{3}^{2}m_{b}n_{a}\theta_{a}^{n_{a}}}{(\theta_{a}^{n_{a}}+{p_{a}^{n_{a}}}^{{SN}^{*}})^{3}}\bigg{\\{}{p_{a}^{n_{a}-2}}^{{SN}^{*}}\theta_{a}^{n_{a}}(n_{a}-1)-(n_{a}+1){p_{a}^{2n_{a}-2}}^{{SN}^{*}}\bigg{\\}}$
($\mathcal{A}_{2}$) $\displaystyle\neq 0,$

where
$\mathcal{X}=\theta_{aa}^{n_{aa}}\theta_{b}^{n_{b}}+{p_{a}^{n_{aa}}}^{{SN}^{*}}\theta_{b}^{n_{b}}+\theta_{aa}^{n_{aa}}{p_{b}^{n_{b}}}^{{SN}^{*}}.$

####  3.1.3 Numerical simulation and discussion corresponding to model (1)

The purpose of this section is to provide validation for the theoretical
results that were covered in Subsection 3.1. We have talked about the system’s
behavior in terms of monostability and bi-stability with fictitious
parameters.

######  Example 1.

$n_{aa}=n_{a}=n_{b}=1$ (case 1) and other parameters are chosen as
$\delta_{a}=0.9,\,\delta_{b}=0.001,\,\gamma_{a}=0.99,\,\gamma_{b}=0.046,\,k_{a}=4.74,\,k_{b}=0.017,\,m_{b}=4.6,\,m_{a}=18.6,\,\theta_{a}=1.8,\,A_{1}=1.1,\,B_{1}=0.01,\,\theta_{aa}=2.6,\,\theta_{b}=(0.1,6.14).$
In this range of $\theta_{b},$ there is unique ($S_{1}^{*}$) steady state in
$\theta_{b}^{1}$, three ($S_{1}^{*},\,S_{2}^{*},\,S_{3}^{*}$) in
$\theta_{b}^{2}$ and unique ($S_{3}^{*}$) in $\theta_{b}^{3},$ such that
$\theta_{b}=\theta_{b}^{1}\cup\theta_{b}^{2}\cup\theta_{b}^{2}$ where
$\theta_{b}^{1}=(0.1,{\theta_{b}^{SN}}^{1}),\,\theta_{b}^{2}=[{\theta_{b}^{SN}}^{1},{\theta_{b}^{SN}}^{2}]$
and $\theta_{b}^{2}=({\theta_{b}^{SN}}^{2},6.14).$ The steady states
$S_{1}^{*},\,S_{3}^{*}$ are stable and ($S_{2}^{*}$) is unstable. The
corresponding bifurcation diagram is shown in Fig. 4 (a).
To ensure the stability of $S_{1}^{*},$ we chosen
$\theta_{b}=1.364\in\theta_{b}^{1}$ and the coefficients of characteristic
equation at $S_{1}^{*}=(1.2963,21.0814,6.8273,358.39)$ are calculated as
$\epsilon_{1}=1.937>0,\,\epsilon_{2}=0.854>0,\,\epsilon_{3}=0.0360>0,\,\epsilon_{4}=0.00003>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.0582>0.$
Therefore, from Theorem 1, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
parameters, shown in Fig. 4 (b). We can see that all trajectories starting
from any initial points are converging to $S_{1}^{*}.$
Further, we chosen $\theta_{b}=3.051\in\theta_{b}^{2}$ and the coefficients of
characteristic equation at $S_{1}^{*}=(1.9213,15.3197,10.1187,260.4357)$ are
calculated as
$\epsilon_{1}=1.9370>0,\,\epsilon_{2}=0.6204>0,\,\epsilon_{3}=0.0251>0,\,\epsilon_{4}=0.0000108>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.0295>0.$
The coefficients of characteristic equation at
$S_{2}^{*}=(4.5514,7.2020,23.9709,122.4344)$ are calculated as
$\epsilon_{1}=1.937>0,\,\epsilon_{2}=0.4297>0,\,\epsilon_{3}=0.0161>0,\,\epsilon_{4}=-0.0000066<0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.0132>0.$
The coefficients of characteristic equation at
$S_{3}^{*}=(11.9929,2.9882,63.1625,50.7998)$ are calculated as
$\epsilon_{1}=1.937>0,\,\epsilon_{2}=0.6397>0,\,\epsilon_{3}=0.0260>0,\,\epsilon_{4}=0.000012>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.0315>0.$
Therefore, from Theorem 1, we conclude that $S_{1}^{*},\,S_{3}^{*}$ is locally
asymptotically stable and $S_{2}^{*}$ is unstable. We draw the phase portrait
for the considered parameters, shown in Fig. 4 (c). We can see that all
trajectories are converging to $S_{1}^{*},\,S_{3}^{*},$ and nearby
trajectories of $S_{2}^{*}$ are moving away from $S_{2}^{*}.$ This shows the
bistable behaviour of the system.
We chosen $\theta_{b}=5.795\in\theta_{b}^{3},$ the coefficients of
characteristic equation at $S_{3}^{*}=(16.3872,$
$2.2604,86.3061,38.4265)$ are calculated as
$\epsilon_{1}=1.937>0,\,\epsilon_{2}=0.8246>0,\,\epsilon_{3}=0.0347>0,\,\epsilon_{4}=0.00002835>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.0541>0.$
Therefore, from Theorem 1, we conclude that $S_{3}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the chosen parameters,
shown in Fig. 4 (d). We can see that all trajectories starting from any
initial points are converging to $S_{3}^{*}.$
For this case, again we choose the parameters as
$\delta_{a}=0.7,\,\delta_{b}=0.1,\,\gamma_{a}=0.67,\,\gamma_{b}=0.265,\,k_{a}=0.9,\,k_{b}=0.06,\,m_{b}=0.9,\,m_{a}=9,\,\theta_{a}=0.86,\,A_{1}=0.8,\,B_{1}=0.9,\,\theta_{aa}=0.86,\,\theta_{b}=(0.00000001,0.4).$
For this set of parameters, we obtain a unique stable steady state as
discussed in Case 1. The corresponding bifurcation plot is shown in Fig. 4
(e). To ensure the stability of $S_{1}^{*},$ we choose $\theta_{b}=0.1$ from
the range of $\theta_{b}.$ The coefficients of characteristic equation at
$S_{1}^{*}$ are calculated as
$\epsilon_{1}=1.735>0,\,\epsilon_{2}=0.738>0,\,\epsilon_{3}=0.11349>0,\,\epsilon_{4}=0.0049>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.1178>0.$
Therefore, from Theorem 1, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
$\theta_{b}$, shown in Fig. 4 (f). We can see that all trajectories starting
from any initial points are converging to
$S_{1}^{*}=(3.61984,3.92592,4.65408,2.35555).$
In this example, we also discuss the case of saddle-node bifurcation. We have
seen that the system (1) has two steady states $S_{2}^{*},\,S_{3}^{*}$ in
$\theta_{b}^{2}$ such that as the value of $\theta_{b}$ decreases, the two
steady states collide at $\theta_{b}^{{{SN}^{1}}}=2.714726411506651$ and
denoted as $S_{{SN}^{1}}^{*}=(8.3509,4.1491,43.9812,70.5353).$ The Jacobian
matrix $J=Df(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})$ has a simple
eigenvalue $\lambda=0$ and $v=(-0.1066,0.0482,-0.5612,0.8194)^{T},$
$w=(-0.0156,0.3466,-0.0033,0.9379)^{T}$, are the eigenvectors of $J$ and
$J^{T},$ respectively. Both the tranversality conditions
$w^{T}f_{\theta_{b}}(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})=-0.0244\neq 0$
and
$w^{T}[D^{2}f(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})(v,v)]=0.0000124\neq 0$
are satisfied. Hence, the system (1) experiences saddle-node bifurcation at
$\theta_{b}=\theta_{b}^{{{SN}^{1}}}.$ Also, the system (1) has two steady
states $S_{1}^{*},\,S_{2}^{*}$ in $\theta_{b}^{2}$ such that as the value of
$\theta_{b}$ increases, the two steady states collide at
$\theta_{b}^{{{SN}^{2}}}=3.345369897729079$ and denoted as
$S_{{SN}^{2}}^{*}=(2.7005,11.4515,14.2227,194.6749).$ The Jacobian matrix
$J=Df(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})$ has a simple eigenvalue
$\lambda=0$ and
$v=(-0.0158,0.0585,-0.0835,0.9947)^{T},\,w=(-0.1278,0.3437,-0.0267,0.93)^{T}$,
are the eigenvectors of $J$ and $J^{T},$ respectively. Both the tranversality
conditions
$w^{T}f_{\theta_{b}}(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})=-0.0541\neq 0$
and
$w^{T}[D^{2}f(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})(v,v)]=-0.000006817\neq
0$ are satisfied. Hence, the system (1) experiences saddle-node bifurcation at
$\theta_{b}=\theta_{b}^{{{SN}^{2}}}.$

Figure 4: For $n_{a}=1,\,n_{b}=1,\,n_{aa}=1$; (a) Bifurcation plot exhibiting
hysteresis effect with three steady states (b) Phase portrait at
$\theta_{b}=1.364$ (c) Phase portrait at $\theta_{b}=3.051$ depicting
bistability (d) Phase portrait at $\theta_{b}=5.795$ (e) Bifurcation plot with
unique stable steady state (f) Phase portrait at $\theta_{b}=0.1.$

######  Example 2.

$n_{a}=n_{b}=1,\,n_{aa}=2$ (case 2), and the other parameters are chosen as:
$\delta_{a}=0.9,\,\delta_{b}=0.04,\,\gamma_{a}=0.604,\,\gamma_{b}=0.3,\,k_{a}=0.525,\,k_{b}=0.9,\,m_{b}=0.01,\,m_{a}=15,\,\theta_{a}=0.4,\,A_{1}=0.65,\,B_{1}=0.785,\,\theta_{aa}=0.7,\,\theta_{b}=(0.1,0.99).$
In this range of $\theta_{b},$ there is unique ($S_{1}^{*}$) steady state in
$\theta_{b}^{1}$, three ($S_{1}^{*},\,S_{2}^{*},\,S_{3}^{*}$) in
$\theta_{b}^{2}$ and unique ($S_{3}^{*}$) in $\theta_{b}^{3},$ such that
$\theta_{b}=\theta_{b}^{1}\cup\theta_{b}^{2}\cup\theta_{b}^{2}$ where
$\theta_{b}^{1}=(0.1,{\theta_{b}^{SN}}^{1}),\,\theta_{b}^{2}=[{\theta_{b}^{SN}}^{1},{\theta_{b}^{SN}}^{2}]$
and $\theta_{b}^{2}=({\theta_{b}^{SN}}^{2},0.99).$ The steady states
$S_{1}^{*},\,S_{3}^{*}$ are stable and ($S_{2}^{*}$) is unstable. The
corresponding bifurcation diagram is shown in Fig. 5 (a).
To ensure the stability of $S_{1}^{*},$ we chosen
$\theta_{b}=0.3\in\theta_{b}^{1}$ and the coefficients of characteristic
equation at $S_{1}^{*}=(1.2011,2.6288,0.7006,59.1476)$ are calculated as
$\epsilon_{1}=1.8440>0,\,\epsilon_{2}=0.9545>0,\,\epsilon_{3}=0.1646>0,\,\epsilon_{4}=0.0052>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.2450>0.$
Therefore, from Theorem 1, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
parameters, shown in Fig. 5 (b). We can see that all trajectories starting
from any initial points are converging to $S_{1}^{*}.$
Further, we chosen $\theta_{b}=0.7\in\theta_{b}^{2}$ and the coefficients of
characteristic equation at $S_{1}^{*}=(1.5546,2.6269,0.9069,59.1046)$ are
calculated as
$\epsilon_{1}=1.8440>0,\,\epsilon_{2}=0.7388>0,\,\epsilon_{3}=0.0913>0,\,\epsilon_{4}=0.0026>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.1073>0.$
The coefficients of characteristic equation at
$S_{2}^{*}=(4.2128,2.6213,2.4575,58.98)$ are calculated as
$\epsilon_{1}=1.8440>0,\,\epsilon_{2}=0.3597>0,\,\epsilon_{3}=-0.0376<0,\,\epsilon_{4}=-0.0020<0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=-0.0196<0.$
The coefficients of characteristic equation at
$S_{3}^{*}=(20.1513,2.6178,11.7549,58.8997)$ are calculated as
$\epsilon_{1}=1.844>0,\,\epsilon_{2}=0.8283>0,\,\epsilon_{3}=0.1217>0,\,\epsilon_{4}=0.0037>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.1587>0.$
Therefore, from Theorem 1, we conclude that $S_{1}^{*},\,S_{3}^{*}$ is locally
asymptotically stable and $S_{2}^{*}$ is unstable. We draw the phase portrait
for the considered parameters, shown in Fig. 5 (c). We can see that all
trajectories are converging to $S_{1}^{*},\,S_{3}^{*},$ and nearby
trajectories of $S_{2}^{*}$ are moving away from $S_{2}^{*}.$ This shows the
bistable behavior of the system.
We chosen $\theta_{b}=0.9\in\theta_{b}^{3},$ the coefficients of
characteristic equation at $S_{3}^{*}=(21.7247,2.6177,$
$12.6727,58.8979)$ are calculated as
$\epsilon_{1}=1.8440>0,\,\epsilon_{2}=0.8928>0,\,\epsilon_{3}=0.1437>0,\,\epsilon_{4}=0.0044>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.2008>0.$
Therefore, from Theorem 1, we conclude that $S_{3}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the chosen parameters,
shown in Fig. 5 (d). We can see that all trajectories starting from any
initial points are converging to $S_{3}^{*}.$
In this example, we also discuss the case of saddle-node bifurcation. We have
seen that the system (1) has two steady states $S_{2}^{*},\,S_{3}^{*}$ in
$\theta_{b}^{2}$ such that as the value of $\theta_{b}$ decreases, the two
steady states collide at $\theta_{b}^{{{SN}^{1}}}=0.466081314925648$ and
denoted as $S_{{SN}^{1}}^{*}=(12.3079,2.6184,7.1796,58.9146).$ The Jacobian
matrix $J=Df(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})$ has a simple
eigenvalue $\lambda=0$ and $v=(0.8638,-1.1694,0.5039,-0.0026)^{T},$
$w=(0.1932,-0.9066,0.2222,-0.3022)^{T}$, are the eigenvectors of $J$ and
$J^{T},$ respectively. Both the tranversality conditions
$w^{T}f_{\theta_{b}}(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})=1.528\neq 0$
and
$w^{T}[D^{2}f(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})(v,v)]=-0.00089\neq 0$
are satisfied. Hence, the system (1) experiences saddle-node bifurcation at
$\theta_{b}=\theta_{b}^{{{SN}^{1}}}.$ Also, the system (1) has two steady
states $S_{1}^{*},\,S_{2}^{*}$ in $\theta_{b}^{2}$ such that as the value of
$\theta_{b}$ increases, the two steady states collide at
$\theta_{b}^{{{SN}^{2}}}=0.845253711971854$ and denoted as
$S_{{SN}^{2}}^{*}=(2.2633,2.6244,1.3202,59.0494).$ The Jacobian matrix
$J=Df(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})$ has a simple eigenvalue
$\lambda=0$ and
$v=(0.8627,-0.0023,0.5032,-0.0510)^{T},\,w=(0.5647,-0.4828,0.6497,-0.1609)^{T}$,
are the eigenvectors of $J$ and $J^{T},$ respectively. Both the tranversality
conditions
$w^{T}f_{\theta_{b}}(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})=0.4497\neq 0$
and $w^{T}[D^{2}f(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})(v,v)]=-0.0081\neq
0$ are satisfied. Hence, the system (1) experiences saddle-node bifurcation at
$\theta_{b}=\theta_{b}^{{{SN}^{2}}}.$
The system (1) exhibits a hysteresis effect in $\theta_{b}^{2}$ where multiple
steady states coexist, as shown in Fig. 4 (a) and Fig 5 (a). The two outer
steady states are stable, while the interior steady state (red) is unstable.
For this case, again we choose
$\gamma_{a}=0.14,\,m_{a}=2,\,A_{1}=0.35,\,\theta_{b}\in(0.1,1.8)$ and other
parameters are the same as above. We obtained a unique steady state that is
stable throughout the range of considered $\theta_{b}.$ The corresponding
bifurcation diagram is shown in Fig. 5 (e). The phase trajectory for
$\theta_{b}=0.3,$ is shown in Fig. 5 (f). All trajectories are converging to
$S_{1}^{*}=(2.91359,2.62302,1.6996,59.0179).$

Figure 5:  For $n_{a}=1,\,n_{b}=1,\,n_{aa}=2$; (a) Bifurcation plot exhibiting
hysteresis effect with three steady states (b) Phase portrait at
$\theta_{b}=0.3$ (c) Phase portrait at $\theta_{b}=0.7,$ depicting bistability
(d) Phase portrait at $\theta_{b}=0.9$ (e) Bifurcation plot with unique stable
steady state (f) Phase portrait at $\theta_{b}=0.3.$

######  Example 3.

$n_{a}=2,\,n_{b}=1,\,n_{aa}=2$ (case 3), and the other parameters are chosen
as:
$\delta_{a}=0.4,\,\delta_{b}=0.001,\,\gamma_{a}=0.67,\,\gamma_{b}=0.046,\,k_{a}=0.7,\,k_{b}=0.017,\,m_{b}=8.375,\,m_{a}=9,\,\theta_{a}=0.6,\,A_{1}=1,\,B_{1}=0.01,\,\theta_{aa}=2.8,\,\theta_{b}\in(0.1,2.9).$
In this range of $\theta_{b},$ there is unique ($S_{3}^{*}$) steady state in
$\theta_{b}^{1}$, three ($S_{1}^{*},\,S_{2}^{*},\,S_{3}^{*}$) in
$\theta_{b}^{2}$ and unique in $\theta_{b}^{3},$ such that
$\theta_{b}=\theta_{b}^{1}\cup\theta_{b}^{2}\cup\theta_{b}^{2}$ where
$\theta_{b}^{1}=(0.1,{\theta_{b}^{SN}}^{1}),\,\theta_{b}^{2}=[{\theta_{b}^{SN}}^{1},{\theta_{b}^{SN}}^{2}]$
and $\theta_{b}^{2}=({\theta_{b}^{SN}}^{2},0.99).$ The steady states
$S_{1}^{*},\,S_{3}^{*}$ are stable and ($S_{2}^{*}$) is unstable. The
corresponding bifurcation diagram is shown in Fig. 6 (a). All other
descriptions are the same as those discussed in Example 2.

######  Example 4.

$n_{a}=2,\,n_{b}=2,\,n_{aa}=1$ (case 4), and the other parameters are chosen
as:
$\delta_{a}=0.9,\,\delta_{b}=0.001,\,\gamma_{a}=1.23,\,\gamma_{b}=0.046,\,k_{a}=5.1,\,k_{b}=0.017,\,m_{b}=10,\,m_{a}=14,\,\theta_{a}=0.6,\,A_{1}=1,\,B_{1}=0.01,\,\theta_{aa}=2.7,\,\theta_{b}\in(0.8,5.8).$
In this range of $\theta_{b},$ there is unique ($S_{3}^{*}$) steady state in
$\theta_{b}^{1}$, three ($S_{1}^{*},\,S_{2}^{*},\,S_{3}^{*}$) in
$\theta_{b}^{2}$ and unique in $\theta_{b}^{3},$ such that
$\theta_{b}=\theta_{b}^{1}\cup\theta_{b}^{2}\cup\theta_{b}^{2}$ where
$\theta_{b}^{1}=(0.1,{\theta_{b}^{SN}}^{1}),\,\theta_{b}^{2}=[{\theta_{b}^{SN}}^{1},{\theta_{b}^{SN}}^{2}]$
and $\theta_{b}^{2}=({\theta_{b}^{SN}}^{2},0.99).$ The steady states
$S_{1}^{*},\,S_{3}^{*}$ are stable and ($S_{2}^{*}$) is unstable. The
corresponding bifurcation diagram is shown in Fig. 6 (b). All other
descriptions are the same as discussed in Example 2.

Figure 6: (a) For $n_{a}=2,\,n_{b}=1,\,n_{aa}=2;$ Bifurcation plot exhibiting
hysteresis effect with three steady states (b) For
$n_{a}=2,\,n_{b}=2,\,n_{aa}=1;$ Bifurcation plot exhibiting hysteresis effect
with three steady states.

Now, we discuss some other cases in which we observe the monostable and
bistable behavior of the system. Instead of showing bifurcation diagrams, here
we show the phase portraits with null clines.
If $n_{a}=2,\,n_{b}=2,\,n_{aa}=2$, and the other parameters are chosen as:
$\delta_{a}=0.2,\,\delta_{b}=0.001,\,\gamma_{a}=0.67,\,\gamma_{b}=.046,\,k_{a}=0.7,\,k_{b}=0.017,\,m_{b}=8.375,\,m_{a}=9,\,\theta_{a}=0.6,\,A_{1}=1,\,B_{1}=0.01,\,\theta_{aa}=2.8,\,\theta_{b}=0.7.$
For this parameter, we obtained three steady states
$S_{1}^{*}=(1.50492,2.5496,5.2672,43.3431),\,S_{2}^{*}=(5.00231,0.430963,17.5081,7.3264)$
and
$S_{3}^{*}=(13.3447,0.247431,46.7065,4.2063)$ which are nodal sink (stable),
saddle point (unstable) and nodal sink, respectively. The corresponding phase
portrait diagram is shown in Fig. 7 (a). Further, for this case, we choose
$\theta_{b}=0.3,$ we obtained unique steady state
$S_{1}^{*}=(1.49473,2.5811,5.23155,43.8786),$ which is locally asymptotically
stable (nodal sink), shown in Fig. 7 (b).

Figure 7: For $n_{a}=n_{b}=n_{aa}=2;$ (a) Phase portrait at $\theta_{b}=0.7,$
depicting bistability (b) Phase portrait at $\theta_{b}=0.3$.

If $n_{a}=3,\,n_{b}=3,\,n_{aa}=3$, and the other parameters are chosen as same
as the case of $n_{a}=2,\,n_{b}=2,\,n_{aa}=2.$ For this parameter, we obtained
three steady states
$S_{1}^{*}=(1.56647,0.455699,5.48266,7.74688),\,S_{2}^{*}=(2.2693,0.295845,7.94254,5.02937)$
and $S_{3}^{*}=(14.6048,0.217686,51.1168,3.70066)$ which are nodal sink
(stable), saddle point (unstable) and nodal sink, respectively. Fig. 8 (a)
shows the corresponding phase portrait diagram. Further, for this case, we
choose $\theta_{b}=0.3,$ we obtained unique steady state
$S_{1}^{*}=(1.49663,0.490592,5.2382,8.34006),$ which is locally asymptotically
stable (nodal sink), shown in Fig. 8 (b).

Figure 8: For $n_{a}=n_{b}=n_{aa}=3;$ (a) Phase portrait at $\theta_{b}=0.7,$
depicting bistability (b) Phase portrait at $\theta_{b}=0.3$.

###  3.2 Conclusion and biological interpretation for this section

In this section, we assume that the two genes $a$ and $b$ compete for binding
to the promoter regions and any of the binding events can elicit gene
expression. This type of biological regulation is modeled using a competitive
OR logic function. Our results show that with any combination and degree of
multimerization (Hill coefficient values), the network will exhibit at most
bistability. We also found that cell state transition can occur through two
mechanisms: (i) hysteresis; in which the cell remembers the past state and
abruptly switches back and forth to another state as the input signal crosses
a certain threshold (Fig.’s 4a, 5a, 6). (ii) biphasic transition; in which we
observe that cells can transit to the next state without a step-like switch
(4c, 5e). In this situation, there is no bifurcating attractor and hence no
cells are ”forced” to transition, instead, gene expression changes in a smooth
continuous way.

##  4 Mathematical Model [AND Logic]

In this section, we model the co-regulation of node $a$ using non-competitive
AND logic [35]. With this formalism, the model system takes the following
form.

$\displaystyle\frac{dr_{a}}{dt}=\frac{m_{a}(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}}{1+(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}(\frac{p_{b}}{\theta_{b}})^{n_{b}}}-\gamma_{a}r_{a}+A_{1},$
(8)
$\displaystyle\frac{dr_{b}}{dt}=\frac{m_{b}}{1+(\frac{p_{a}}{\theta_{a}})^{n_{a}}}-\gamma_{b}r_{b}+B_{1},$
$\displaystyle\frac{dp_{a}}{dt}=k_{a}r_{a}-\delta_{a}p_{a},$
$\displaystyle\frac{dp_{b}}{dt}=k_{b}r_{b}-\delta_{b}p_{b},$

with $r_{a}(0)\geq 0,\,r_{b}(0)\geq 0,\,p_{a}(0)\geq 0,\,p_{b}\geq 0.$ The
parameter description is same as Section 3.

###  4.1 Mathematical Analysis

This subsection again presents the existence and stability analysis of the
steady states for the model system (8). We use the Routh-Hurwitz criterion to
prove the local asymptotic stability of the steady states, while Sotomayor’s
theorem is used to derive the transversality conditions of saddle-node
bifurcation.

We obtained the steady state $S^{*}=(r_{a}^{*},r_{b}^{*},p_{a}^{*},p_{b}^{*})$
by setting the right hand side of the model (1) to zero. From third and fourth
equations of model (1), we obtained
${r_{a}}^{*}=\frac{\delta_{a}{p_{a}}^{*}}{k_{a}}$ and
${r_{b}}^{*}=\frac{\delta_{b}{p_{b}}^{*}}{k_{b}},$ respectively, as
$\frac{dp_{a}}{dt}=0$ and $\frac{dp_{b}}{dt}=0.$ Using $\frac{dr_{b}}{dt}=0$
and then putting the value of $r_{b},$ we get
${p_{b}}^{*}=\frac{1}{\gamma_{b}}\bigg{\\{}B_{1}^{\prime}+\frac{m_{2}}{\theta_{a}^{n_{a}}+{{p_{a}}^{*}}^{n_{a}}}\bigg{\\}}.$
Finally, using the value of $p_{b}$ and $r_{a}$ in $\frac{dr_{a}}{dt}=0,$ we
get a polynomial in $p_{a}$ as

$\gamma_{b}^{n_{b}}(\theta_{a}^{n_{a}}+p_{a}^{n_{a}})^{n_{b}}\bigg{(}m_{1}p_{a}^{n_{aa}}+\theta_{b}^{n_{b}}\theta_{aa}^{n_{aa}}(A_{1}^{\prime}-\gamma_{a}p_{a})\bigg{)}+(A_{1}^{\prime}-\gamma_{a}p_{a})p_{a}^{n_{aa}}\bigg{(}m_{2}+B_{1}^{\prime}(\theta_{a}^{n_{a}}+p_{a}^{n_{a}})\bigg{)}^{n_{b}}=0,$
(9)

where
$A_{1}^{\prime}=\frac{A_{1}k_{a}}{\delta_{a}},\,B_{1}^{\prime}=\frac{B_{1}k_{b}}{\delta_{b}},\,m_{1}=\frac{m_{a}k_{a}\theta_{b}^{n_{b}}}{\delta_{a}},$
and $m_{2}=\frac{m_{b}k_{b}\theta_{a}^{n_{a}}}{\delta_{b}}.$ Here, $p_{a}^{*}$
is a positive real root of (9) which we discuss in different cases explicitly.
Case 1: When $n_{aa}=n_{a}=n_{b}=1,$ then from (9), we get a cubic polynomial
in $p_{a}$

$f(p_{a})=Ap_{a}^{3}+Bp_{a}^{2}+Cp_{a}+D,$ (10)

where
$A=-\gamma_{a}B_{1}^{\prime},\,B=\gamma_{b}(m_{1}-\gamma_{a}\theta_{aa}\theta_{b})-\gamma_{a}(m_{2}+B_{1}^{\prime}\theta_{a})+A_{1}^{\prime}B_{1}^{\prime}\\\
C=\gamma_{b}A_{1}^{\prime}\theta_{b}\theta_{aa}+\gamma_{b}\theta_{a}(m_{1}-\gamma_{a}\theta_{aa}\theta_{b})+A_{1}^{\prime}(m_{2}+B_{1}^{\prime}\theta_{a})$
and $D=\gamma_{b}\theta_{a}A_{1}^{\prime}\theta_{b}\theta_{aa}.$ The nature of
the curve $f$ can be observe from (10) as

  1. (i)

$f(0)>0$ because $D>0,$

  2. (ii)

$f\to-\infty\,(\infty)$ as $p_{a}\to\infty\,(-\infty)$ because $A<0,$

  3. (iii)

The number of positive real roots of equation (10) can be seen using
Descarte’s rule of sign, mentioned in below Table (3).

Coefficients |  Number of possible positive real roots of $f(p_{a})=0$
---|---
A | B | C | D
- | - | - | + | 1
- | + | - | + | 3,1
- | + | + | + | 1
- | - | + | + | 1
Table 3: Number of possible positive real roots of $f(p_{a})=0.$

Therefore, $f(p_{a})=0$ will have either unique or three positive real roots.
Now, to discard the existence of three positive real roots, we use the
nullcline method. From the system (8), we have

$\displaystyle
F_{1}(p_{a},p_{b}):=\frac{m_{a}k_{a}\theta_{b}p_{a}}{\delta_{a}(\theta_{aa}\theta_{b}+p_{a}p_{b})}+A_{1}^{\prime}-\gamma_{a}p_{a}$
(11) $\displaystyle
F_{2}(p_{a},p_{b}):=\frac{m_{b}k_{b}\theta_{a}}{\delta_{b}(\theta_{a}+p_{a})}+B_{1}^{\prime}-\gamma_{b}p_{b},$

First to observe the nature of curve $F_{2}(p_{a},p_{b})=0,$ we reduce it in
the form of

$p_{b}=g_{1}(p_{a}):=\frac{1}{k_{2}}\Big{\\{}\frac{m_{b}\theta_{a}}{(\theta_{a}+p_{a})}+B_{1}\Big{\\}}.$
($\mathscr{A}_{1}$)

The following conclusion can be made from ($\mathscr{A}_{1}$)

  1. $\mathscr{A}_{1}$(i)

$g_{1}(0)=\frac{1}{k_{2}}(m_{b}+B_{1})>\frac{B_{1}}{k_{2}}>0$

  2. $\mathscr{A}_{1}$(ii)

$g_{1}(p_{a})\to\frac{B_{1}}{k_{2}}$ as $p_{a}\to\pm\infty$

  3. $\mathscr{A}_{1}$(iii)

$g^{\prime}_{1}(p_{a})<0$, i.e, $g_{1}(p_{a})$ is a decreasing function of
$p_{a}.$

The graph of $F_{2}(p_{a},p_{b})=0$ (or $p_{b}=g_{1}(p_{a})$) is shown in Fig.
9.
Now, to observe the nature of curve $F_{1}(p_{a},p_{b})=0,$ we reduce it in
the form of

$p_{b}=g_{2}(p_{a}):=\frac{\theta_{b}(m_{a}p_{a}+(A_{1}-k_{1}p_{a})\theta_{aa})}{p_{a}(k_{1}p_{a}-A_{1})}.$
($\mathscr{B}_{1}$)

Equivalently

$p_{b}=g_{2}(p_{a}):=\frac{p_{a}\big{(}m_{a}-k_{1}\theta_{aa}\big{)}+A_{1}k_{1}}{p_{a}(k_{1}p_{a}-A_{1})}.$
($\mathscr{B}_{2}$)

Also, we find

$g_{2}^{\prime}(p_{a})=\frac{-\theta_{b}\Big{(}k_{1}p_{a}^{2}\big{(}m_{a}-A_{1}\theta_{aa}\big{)}+k_{1}p_{a}A_{1}\theta_{aa}+A_{1}\theta_{aa}\big{(}k_{1}p_{a}-A_{1}\big{)}\Big{)}}{p_{a}^{2}(k_{1}p_{a}-A_{1})^{2}}.$
($\mathscr{C}_{1}$)

Equivalently

$g_{2}^{\prime}(p_{a})=\frac{-\theta_{b}\bigg{(}k_{1}p_{a}\Big{(}p_{a}(m_{a}-k_{1}\theta_{aa})+A_{1}\theta_{aa}\Big{)}+A_{1}\theta_{aa}\big{(}k_{1}p_{a}-A_{1}\big{)}\bigg{)}}{p_{a}^{2}(k_{1}p_{a}-A_{1})^{2}}.$
($\mathscr{C}_{2}$)

The following conclusion can be made from ($\mathscr{B}_{1}$)

  1. $\mathscr{B}$(i)

$g_{2}(p_{a})$ is a rational function, and it has two vertical asymptotes at
$p_{a}=0$ and $p_{a}=\frac{A_{1}}{k_{1}}$ and $p_{b}=0$ (i.e., $p_{a}$-axis)
is a horizontal asymptote.

  2. $\mathscr{B}$(ii)

We note that for $p_{a}\in\big{(}0,\frac{A_{1}}{k_{1}}\big{)}$, we have
$g_{2}(p_{a})<0$ so if positive real roots of $g_{2}(p_{a})=0$ exits then it
always lie on the right-hand side of vertical asymptote
$p_{a}=\frac{A_{1}}{k_{1}}$.

  3. $\mathscr{B}$(iii)

The $p_{a}-$intercept can be find with $g_{2}(p_{a})=0,$ which implies that
$p_{a}=-\frac{A_{1}\theta_{aa}}{(m_{a}-k_{1}\theta_{aa})}.$ The sign of
$(m_{a}-k_{1}\theta_{aa})$ will provide the unique $p_{a}-$intercept of
$g_{2}(p_{a}).$ Now, we will discuss three cases to understand the nature of
curve $p_{b}=g_{2}(p_{a}).$

     * (a)

When $m_{a}-k_{1}\theta_{aa}>0$ and $p_{a}>\frac{A_{1}}{k_{1}}$ then there
will no positive intercept; also $g_{2}(p_{a})>0$ and
$g_{2}^{\prime}(p_{a})<0$ which implies that $g_{2}(p_{a})$ is strictly
decreasing. The corresponding diagram is shown in Fig. 10 (a).

     * (b)

When $m_{a}-k_{1}\theta_{aa}<0$ then there will unique positive intercept
$p_{a}=-\frac{A_{1}\theta_{aa}}{(m_{a}-k_{1}\theta_{aa})}.$ 10 (c). When
$\frac{A_{1}}{k_{1}}<p_{a}<-\frac{A_{1}\theta_{aa}}{(m_{a}-k_{1}\theta_{aa})}$
then $g_{2}(p_{a})>0$ and $g_{2}^{\prime}(p_{a})<0$ which implies that
$g_{2}(p_{a})$ is strictly decreasing. The corresponding diagram is shown in
Fig. 10 (c). When $p_{a}>-\frac{A_{1}\theta_{aa}}{(m_{a}-k_{1}\theta_{aa})}$
then $g_{2}(p_{a})<0,$ and never cross the $p_{a}$-axis.

All possibilities of the existence of a steady state are shown in Fig. 10 (b)
and Fig. 10 (d) corresponding to 10 (a) and Fig. 10 (c), respectively.
Therefore, from the above discussion, we discard the existence of three
positive steady states.

Figure 9: For $n_{a}=n_{b}=n_{aa}=1;$ Plot of $g_{1}(p_{a})$ as a decreasing
function of $p_{a}.$

Figure 10: For $n_{a}=n_{b}=n_{aa}=1;$ (a) The graph depicting no positive
$p_{a}$-intercept of $g_{2}(p_{a})$ (b) Existence of unique positive steady
state corresponding to Fig. 10 (a). (c) The graph depicting unique positive
$p_{a}$-intercept of $g_{2}(p_{a})$ (d) Existence of unique positive steady
state corresponding to Fig. 10 (c).

Case 2: When $n_{a}=n_{b}=1,\,n_{aa}=2,$ then from (8), we get a quartic
polynomial in $p_{a}$

$f(p_{a})=Ap_{a}^{4}+Bp_{a}^{3}+Cp_{a}^{2}+Dp_{a}+E,$ (12)

where
$A=B_{1}^{\prime}\,\gamma_{a},\,B=-A_{1}^{\prime}\,B_{1}^{\prime}-\gamma_{b}\,m_{1}+\gamma_{a}\,\left(m_{2}+B_{1}^{\prime}\,\theta_{a}\right),\,C=\gamma_{a}\,\gamma_{b}\,\theta_{b}\,{\theta_{aa}}^{2}-A_{1}^{\prime}\,\left(m_{2}+B_{1}^{\prime}\,\theta_{a}\right)\\\
-\gamma_{b}\,m_{1}\,\theta_{a},\,D=-\gamma_{b}\,\theta_{b}\,{\theta_{aa}}^{2}\,\left(A_{1}^{\prime}-\gamma_{a}\,\theta_{a}\right),$
and
$E=-A_{1}^{\prime}\,\gamma_{b}\,\theta_{a}\,\theta_{b}\,{\theta_{aa}}^{2}.$
The number of positive real roots of equation (12) can be seen using
Descarte’s rule of sign (mentioned in below Table (2)); therefore,
$f(p_{a})=0$ will have either unique or three positive real roots.

####  4.1.1 Local stability of steady state

######  Theorem 2.

The steady state $S_{1}^{*}$ is locally asymptotically stable if
$\epsilon_{1}>0,\,\epsilon_{3}>0,\,\epsilon_{4}>0,$ and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}>0,$
where $\epsilon_{i},\,i=1,..,4$ are provided in the proof.

###### Proof.

At $S_{1}^{*},$ the Jacobian matrix of the model system (1) is given by

$J_{S_{1}^{*}}=\left[\begin{array}[]{cccc}-\gamma_{a}&0&X&Y\\\
0&-\gamma_{b}&Z&0\\\ k_{a}&0&-\delta_{a}&0\\\
0&-k_{b}&0&-\delta_{b}\end{array}\right].$

where
$X=\frac{\mathrm{m_{a}}\,n_{aa}\,{\mathrm{p_{a}}}^{n_{aa}-1}\,{\theta_{b}}^{2\,n_{b}}\,{\theta_{aa}}^{n_{aa}}}{{\left({\mathrm{p_{a}}}^{n_{aa}}\,{\mathrm{p_{b}}}^{n_{b}}+{\theta_{b}}^{n_{b}}\,{\theta_{aa}}^{n_{aa}}\right)}^{2}},\,Y=-\frac{\mathrm{m_{a}}\,n_{b}\,{\mathrm{p_{a}}}^{2\,n_{aa}}\,{\mathrm{p_{b}}}^{n_{b}-1}\,{\theta_{b}}^{n_{b}}}{{\left({\mathrm{p_{a}}}^{n_{aa}}\,{\mathrm{p_{b}}}^{n_{b}}+{\theta_{b}}^{n_{b}}\,{\theta_{aa}}^{n_{aa}}\right)}^{2}},$
and
$Z=-\frac{\mathrm{m}_{b}\,n_{a}\,{\mathrm{p_{a}}}^{n_{a}-1}\,{\theta_{a}}^{n_{a}}}{{\left({\mathrm{p_{a}}}^{n_{a}}+{\theta_{a}}^{n_{a}}\right)}^{2}}.$
The characteristic equation of $J_{S_{1}^{*}}$ is given as

$\lambda^{4}+\epsilon_{1}\lambda^{3}+\epsilon_{2}\lambda^{2}+\epsilon_{3}\lambda+\epsilon_{4}=0.$
(13)

Here,
$\epsilon_{1}=\gamma_{a}+\gamma_{b}+\delta_{a}+\delta_{b},\epsilon_{2}=\gamma_{a}\gamma_{b}+(\gamma_{a}+\gamma_{b})(\delta_{a}+\delta_{b})+\delta_{a}\delta_{b}-k_{a}X,\,\epsilon_{3}=\delta_{a}\delta_{b}(\gamma_{a}+\gamma_{b})+(\delta_{a}+\delta_{b})(\gamma_{a}\gamma_{b})-k_{a}(\delta_{b}+\gamma_{b})X,$
and $\epsilon_{4}=\gamma_{a}\gamma_{b}\delta_{a}\delta_{b}-k_{a}k_{b}YZ-
k_{a}\gamma_{b}\delta_{b}X.$ By the Routh-Hurwitz criteria, all the roots of
(13) have negative real parts provided
$\epsilon_{1}>0,\,\epsilon_{3}>0,\,\epsilon_{4}>0,$ and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}>0.$
∎

####  4.1.2 Saddle Node Bifurcation

We find the transversality condition the same as discussed in Subsection
3.1.2. Define $f=(f_{1},f_{2},f_{3},f_{4})^{T},$ where $f_{1}$ are defined as

$f_{1}=\frac{m_{a}(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}}{1+(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}(\frac{p_{b}}{\theta_{b}})^{n_{b}}}-\gamma_{a}r_{a}+A_{1},$

and $f_{2},\,f_{3},$ and $f_{4}$ expression is same as Subsection 3.1.2. Then
the model system (8) experiences a saddle-node bifurcation at
$(S_{SN}^{*},\theta_{b}^{SN})$, if the following transversality conditions
[36] are satisfied:

$w^{T}f_{\theta_{b}}(S_{SN}^{*},\theta_{b}^{SN})=w_{1}\frac{m_{a}n_{b}{\theta_{b}^{n_{b}-1}}^{SN}{p_{b}^{n_{b}}}^{{SN}^{*}}{p_{a}^{2n_{aa}}}^{{SN}^{*}}}{(\mathcal{X})^{2}}\neq
0.$ ($\mathcal{T}_{11}$)

and

$\displaystyle w^{T}[D^{2}f(S_{SN}^{*},\theta_{b}^{SN})(v,v)]=$
$\displaystyle\frac{w_{1}}{\mathcal{X}^{3}}\Biggl{\\{}\mathcal{O}\bigg{[}v_{3}^{2}\Big{(}\theta_{aa}^{n_{aa}}{\theta_{b}^{n_{b}}}^{SN}(n_{aa}-1){p_{a}^{(n_{aa}-2)}}^{{SN}^{*}}-(n_{aa}+1){p_{b}^{n_{b}}}^{{SN}^{*}}$
$\displaystyle{p_{a}^{(2n_{aa}-2)}}^{{SN}^{*}}\Big{)}-4v_{3}v_{4}n_{b}{p_{a}^{(2n_{aa}-1)}}^{{SN}^{*}}{p_{b}^{(n_{b}-1)}}^{{SN}^{*}}-v_{4}^{2}\mathcal{O}_{1}$
$\displaystyle\Big{(}\theta_{aa}^{n_{aa}}{\theta_{b}^{n_{b}}}^{SN}(n_{b}-1){p_{b}^{(n_{b}-2)}}^{{SN}^{*}}-(n_{b}+1){p_{a}^{n_{aa}}}^{{SN}^{*}}{p_{b}^{(2n_{b}-2)}}^{{SN}^{*}}\Big{)}\bigg{]}\Biggr{\\}}$
$\displaystyle-\frac{w_{2}v_{3}^{2}m_{b}n_{a}\theta_{a}^{n_{a}}}{(\theta_{a}^{n_{a}}+{p_{a}^{n_{a}}}^{{SN}^{*}})^{3}}\bigg{\\{}{p_{a}^{n_{a}-2}}^{{SN}^{*}}\theta_{a}^{n_{a}}(n_{a}-1)-(n_{a}+1){p_{a}^{2n_{a}-2}}^{{SN}^{*}}\bigg{\\}}$
$\displaystyle\neq 0,$ ($\mathcal{T}_{22}$)

where
$\mathcal{X}=\theta_{aa}^{n_{aa}}\theta_{b}^{n_{b}}+{p_{a}^{n_{aa}}}^{{SN}^{*}}{p_{b}^{n_{b}}}^{{SN}^{*}},\,\mathcal{O}=n_{aa}m_{a}\theta_{b}^{2n_{b}}\theta_{aa}^{n_{aa}},$
and
$\mathcal{O}_{1}=m_{a}n_{b}\theta_{b}^{n_{b}}{p_{aa}^{2n_{aa}}}^{{SN}^{*}}.$

####  4.1.3 Numerical simulation and discussion corresponding to model system
(8)

The purpose of this section is to provide validation for the theoretical
results that were covered in Subsection 4.1. We have talked about the system’s
behavior in terms of monostability and bi-stability with fictitious
parameters.

######  Example 5.

$=n_{a}=n_{b}=n_{aa}=1$ (case 1) and the other parameters are chosen as
$\delta_{a}=0.99,\,\delta_{b}=0.06,\,\gamma_{a}=3.9,\,\gamma_{b}=0.35,\,k_{a}=0.5,\,k_{b}=0.09,\,m_{a}=0.02,\,m_{b}=0.9,\,\theta_{a}=0.0086,\,A_{1}=0.8,\,B_{1}=0.009,\,\theta_{aa}=0.086,\,\theta_{b}=(0.0001,5).$
For this set of parameters, we obtain a unique stable steady state as
discussed in Case 1(and). The corresponding bifurcation plot is shown in Fig.
11 (a). To ensure the stability of $S_{1}^{*},$ we choose $\theta_{b}=2.5$
from the range of $\theta_{b}.$ The coefficients of characteristic equation at
$S_{1}^{*}$ are calculated as
$\epsilon_{1}=5.3>0,\,\epsilon_{2}=5.8008>0,\,\epsilon_{3}=1.6504>0,\,\epsilon_{4}=0.0790>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=45.795>0.$
Therefore, from Theorem 2, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
$\theta_{b}$, shown in Fig. 11 (b). We can see that all trajectories starting
from any initial points are converging to
$S_{1}^{*}=(0.2106,0.2181,0.1064,0.3271).$

Figure 11: For $n_{a}=1,\,n_{b}=1,\,n_{aa}=1$; (a) Bifurcation plot with
unique stable steady state (b) Phase portrait at $\theta_{b}=2.5.$

######  Example 6.

$n_{a}=n_{b}=1,\,n_{aa}=2$ (case 2), and the other parameters are chosen as:
$\delta_{a}=0.8,\,\delta_{b}=0.03,\,\gamma_{a}=0.91,\,\gamma_{b}=0.235,\,k_{a}=0.1,\,k_{b}=0.04,\,m_{b}=0.2,\,m_{a}=6,\,\theta_{a}=0.55,\,A_{1}=5.8,\,B_{1}=0.02,\,\theta_{aa}=1.6,\,\theta_{b}=(0.4,18).$
In this range of $\theta_{b},$ there is unique ($S_{1}^{*}$) steady state in
$\theta_{b}^{1}$, three ($S_{1}^{*},\,S_{2}^{*},\,S_{3}^{*}$) in
$\theta_{b}^{2}$ and unique ($S_{3}^{*}$) in $\theta_{b}^{3},$ such that
$\theta_{b}=\theta_{b}^{1}\cup\theta_{b}^{2}\cup\theta_{b}^{2}$ where
$\theta_{b}^{1}=(0.4,{\theta_{b}^{SN}}^{1}),\,\theta_{b}^{2}=[{\theta_{b}^{SN}}^{1},{\theta_{b}^{SN}}^{2}]$
and $\theta_{b}^{2}=({\theta_{b}^{SN}}^{2},18).$ The steady states
$S_{1}^{*},\,S_{3}^{*}$ are stable and ($S_{2}^{*}$) is unstable. The
corresponding bifurcation diagram is shown in Fig. 12 (a).
To ensure the stability of $S_{1}^{*},$ we chosen
$\theta_{b}=0.6\in\theta_{b}^{1}$ and the coefficients of characteristic
equation at $S_{1}^{*}=(8.4987,0.3754,1.0623,0.5006)$ are calculated as
$\epsilon_{1}=1.975>0,\,\epsilon_{2}=0.922>0,\,\epsilon_{3}=0.1344>0,\,\epsilon_{4}=0.0031>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.2147>0.$
Therefore, from Theorem 2, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
parameters, shown in Fig. 12 (b). We can see that all trajectories starting
from any initial points are converging to $S_{1}^{*}.$
Further, we chosen $\theta_{b}=6\in\theta_{b}^{2}$ and the coefficients of
characteristic equation at $S_{1}^{*}=(10.9693,0.3288,1.3712,0.4383)$ are
calculated as
$\epsilon_{1}=1.975>0,\,\epsilon_{2}=0.6093>0,\,\epsilon_{3}=0.0516>0,\,\epsilon_{4}=0.000992>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.0555>0.$
The coefficients of characteristic equation at
$S_{2}^{*}=(17.4281,0.25666,2.17851,0.342213)$ are calculated as
$\epsilon_{1}=1.975>0,\,\epsilon_{2}=0.3530>0,\,\epsilon_{3}=-0.01634<0,\,\epsilon_{4}=-0.000922<0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=-0.0081<0.$
The coefficients of characteristic equation at
$S_{3}^{*}=(282.698,0.0981496,35.3372,0.130866)$ are calculated as
$\epsilon_{1}=1.975>0,\,\epsilon_{2}=1.0659>0,\,\epsilon_{3}=0.1726>0,\,\epsilon_{4}=0.0036>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.319201502799312>0.$
Therefore, from Theorem 2, we conclude that $S_{1}^{*},\,S_{3}^{*}$ is locally
asymptotically stable and $S_{2}^{*}$ is unstable. We draw the phase portrait
for the considered parameters, shown in Fig. 12 (c). We can see that all
trajectories are converging to $S_{1}^{*},\,S_{3}^{*},$ and nearby
trajectories of $S_{2}^{*}$ are moving away from $S_{2}^{*}.$ This shows the
bistable behaviour of the system.
We chosen $\theta_{b}=17\in\theta_{b}^{3},$ the coefficients of characteristic
equation at $S_{3}^{*}=(924.286,\\\ 0.0891386,115.536,0.118851)$ are
calculated as
$\epsilon_{1}=1.975>0,\,\epsilon_{2}=1.1496>0,\,\epsilon_{3}=0.1947>0,\,\epsilon_{4}=0.0046>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.3861>0.$
Therefore, from Theorem 2, we conclude that $S_{3}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the chosen parameters,
shown in Fig. 12 (d). We can see that all trajectories starting from any
initial points are converging to $S_{3}^{*}.$
In this example, we also discuss the case of saddle-node bifurcation. We have
seen that the system (8) has two steady states $S_{2}^{*},\,S_{3}^{*}$ in
$\theta_{b}^{2}$ such that as the value of $\theta_{b}$ decreases, the two
steady states collide at $\theta_{b}^{{{SN}^{1}}}=2.389134923929666$ and
denoted as $S_{{SN}^{1}}^{*}=(45.7614,0.159759,5.72018,0.21301).$ The Jacobian
matrix $J=Df(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})$ has a simple
eigenvalue $\lambda=0$ and $v=(-0.9923,0.0015,-0.1240,0.002)^{T},$
$w=(-0.00033,0.1678,-0.0030,0.9858)^{T}$, are the eigenvectors of $J$ and
$J^{T},$ respectively. Both the tranversality conditions
$w^{T}f_{\theta_{b}}(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})=-0.0026\neq 0$
and
$w^{T}[D^{2}f(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})(v,v)]=0.00000345\neq
0$ are satisfied. Hence, the system (1) experiences saddle-node bifurcation at
$\theta_{b}=\theta_{b}^{{{SN}^{1}}}.$ Also, the system (8) has two steady
states $S_{1}^{*},\,S_{2}^{*}$ in $\theta_{b}^{2}$ such that as the value of
$\theta_{b}$ increases, the two steady states collide at
$\theta_{b}^{{{SN}^{2}}}=15.667663989915999$ and denoted as
$S_{{SN}^{2}}^{*}=(12.9963,0.300364,1.62454,0.400485).$ The Jacobian matrix
$J=Df(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})$ has a simple eigenvalue
$\lambda=0$ and $v=(-0.9921,\\\
0.0123,-0.124,0.0164)^{T},\,w=(-0.0627,0.1374,-0.5705,0.8073)^{T}$, are the
eigenvectors of $J$ and $J^{T},$ respectively. Both the tranversality
conditions
$w^{T}f_{\theta_{b}}(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})=-0.00062\\\
\neq 0$ and
$w^{T}[D^{2}f(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})(v,v)]=-0.0040\neq 0$
are satisfied. Hence, the system (8) experiences saddle-node bifurcation at
$\theta_{b}=\theta_{b}^{{{SN}^{2}}}.$
The system (8) exhibits a hysteresis effect in $\theta_{b}^{2}$ where multiple
steady states coexist, as shown in Fig 12 (a). The two outer steady states are
stable, while the interior steady state (red) is unstable.
Again, in this case, the parameters are chosen as
$\delta_{a}=0.1,\,\delta_{b}=0.03,\,\gamma_{a}=0.99,\,\gamma_{b}=0.005,\,k_{a}=0.525,\,k_{b}=0.9,\,m_{b}=0.01,\,m_{a}=5,\,\theta_{a}=0.4,\,A_{1}=0.25,\,B_{1}=0.00285,\,\theta_{aa}=20.7,\,\theta_{b}=(0.001,1.2).$
For this set of parameters, we obtain a unique stable steady state as
discussed in Case 2 (AND). The corresponding bifurcation plot is shown in Fig.
12 (e). To ensure the stability of $S_{1}^{*},$ we choose $\theta_{b}=0.9$
from the range of $\theta_{b}.$ The coefficients of characteristic equation at
$S_{1}^{*}$ are calculated as
$\epsilon_{1}=1.1250>0,\,\epsilon_{2}=0.1243>0,\,\epsilon_{3}=0.0032>0,\,\epsilon_{4}=0.00001284>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.0004171>0.$
Therefore, from Theorem 2, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
$\theta_{b}$, shown in Fig. 12 (f). We can see that all trajectories starting
from any initial points are converging to
$S_{1}^{*}=(0.273443,1.00583,1.43558,30.1749).$

Figure 12: For $n_{a}=1,\,n_{b}=1,\,n_{aa}=2$; (a) Bifurcation plot exhibiting
hysteresis effect with three steady states (b) Phase portrait at
$\theta_{b}=0.6$ (c) Phase portrait at $\theta_{b}=6$ depicting bistablility
(d) Phase portrait at $\theta_{b}=17$ (e) Bifurcation plot with the unique
stable steady state (f) Phase portrait at $\theta_{b}=0.9.$

###  4.2 Conclusion and biological interpretation for this section

In this section, we assume that the two genes $a$ and $b$ don’t compete for
binding to the promoter regions. However, only both binding events can elicit
gene expression. This type of biological regulation is modeled using a non-
competitive AND logic function. Our results show that, unlike competitive
binding with OR logic, AND logic can’t produce bistability with monomeric
regulations (Hill coefficients of degree 1). We show that bistability requires
dimeric (Hill coefficients of degree 2) autoloop. We also show that only
through ”forced” cell state transition (hysteresis) and not the smooth state
swap (biphasic transition) is possible with this logic function for all
combinations of dimerizations.

##  5 Mathematical Model [Type-II $\mathcal{OR}$ Logic]

In this section, we model the co-regulation of node $a$ using non-competitive
OR logic [35]. With this formalism, the model system takes the following form.

$\displaystyle\frac{dr_{a}}{dt}=m_{a}\biggl{(}\frac{(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}}{1+(\frac{p_{a}}{\theta_{aa}})^{n_{aa}}}+\frac{1}{1+(\frac{p_{b}}{\theta_{b}})^{n_{b}}}\biggr{)}-\gamma_{a}r_{a}+A_{1},$
(14)
$\displaystyle\frac{dr_{b}}{dt}=\frac{m_{b}}{1+(\frac{p_{a}}{\theta_{a}})^{n_{a}}}-\gamma_{b}r_{b}+B_{1},$
$\displaystyle\frac{dp_{a}}{dt}=k_{a}r_{a}-\delta_{a}p_{a},$
$\displaystyle\frac{dp_{b}}{dt}=k_{b}r_{b}-\delta_{b}p_{b},$

with $r_{a}(0)\geq 0,\,r_{b}(0)\geq 0,\,p_{a}(0)\geq 0,\,p_{b}\geq 0.$
Parameter description is same as Section 3.

###  5.1 Mathematical Analysis

We get the steady state $S^{*}=(r_{a}^{*},r_{b}^{*},p_{a}^{*},p_{b}^{*})$ by
setting the right hand side of the model (14) to zero. From third and fourth
equations of model (14), we obtained
${r_{a}}^{*}=\frac{\delta_{a}{p_{a}}^{*}}{k_{a}}$ and
${r_{b}}^{*}=\frac{\delta_{b}{p_{b}}^{*}}{k_{b}},$ respectively, as
$\frac{dp_{a}}{dt}=0$ and $\frac{dp_{b}}{dt}=0.$ Using $\frac{dr_{b}}{dt}=0$
and then putting the value of $r_{b},$ we get
${p_{b}}^{*}=\frac{1}{\gamma_{b}}\bigg{\\{}B_{1}^{\prime}+\frac{m_{2}}{\theta_{a}^{n_{a}}+{{p_{a}}^{*}}^{n_{a}}}\bigg{\\}}.$
Finally, using the value of $p_{b}$ and $r_{a}$ in $\frac{dr_{a}}{dt}=0,$ we
get a polynomial in $p_{a}$ as

$\displaystyle\bigg{(}m_{2}+$ $\displaystyle
B_{1}^{\prime}\Big{(}\theta_{a}^{n_{a}}+p_{a}^{n_{a}}\Big{)}\bigg{)}^{n_{b}}\biggl{\\{}-\gamma_{a}p_{a}^{n_{aa+1}}+p_{a}^{n_{aa}}\big{(}m_{3}+A_{1}^{\prime}\big{)}-p_{a}\gamma_{a}\theta_{aa}^{n_{aa}}+A_{1}^{\prime}\theta_{aa}^{n_{aa}}\biggr{\\}}+\theta_{b}^{n_{b}}\gamma_{b}^{n_{b}}$
(15)
$\displaystyle(\theta_{a}^{n_{a}}+p_{a}^{n_{a}})^{n_{b}}\biggl{\\{}-\gamma_{a}p_{a}^{n_{aa+1}}+p_{a}^{n_{aa}}\big{(}2m_{3}+A_{1}^{\prime}\big{)}-p_{a}\gamma_{a}\theta_{aa}^{n_{aa}}+(m_{3}+A_{1}^{\prime})\theta_{aa}^{n_{aa}}\biggr{\\}}=0,$

where
$A_{1}^{\prime}=\frac{A_{1}k_{a}}{\delta_{a}},\,B_{1}^{\prime}=\frac{B_{1}k_{b}}{\delta_{b}},\,m_{2}=\frac{m_{b}k_{b}\theta_{a}^{n_{a}}}{\delta_{b}},$
and $m_{3}=\frac{m_{a}k_{a}}{\delta_{a}}.$ Here, $p_{a}^{*}$ is a positive
real root of (15) which we discuss in different cases explicitly.
Case 1: When $n_{a}=n_{b}=n_{aa}=1,$ then from (15), we get a cubic polynomial
in $p_{a}$

$f(p_{a})=Ap_{a}^{3}+Bp_{a}^{2}+Cp_{a}+D,$ (16)

where
$A=-\gamma_{a}(B_{1}^{\prime}+\gamma_{b}\theta_{b}),\,B=-\gamma_{a}(B_{1}^{\prime}\theta_{a}+m_{2})+B_{1}^{\prime}(m_{3}+A_{1}^{\prime}-\gamma_{a}\theta_{aa})-\gamma_{a}\gamma_{b}\theta_{a}\theta_{b}+\theta_{b}\gamma_{b}(2m_{3}+A_{1}^{\prime}-\gamma_{a}\theta_{aa}),\,C=A_{1}^{\prime}B_{1}^{\prime}\theta_{aa}+(B_{1}^{\prime}\theta_{a}+m_{2})(m_{3}+A_{1}^{\prime}-\gamma_{a}\theta_{aa})+\theta_{b}\gamma_{b}\theta_{aa}(m_{3}+A_{1}^{\prime})+\theta_{a}\gamma_{b}\theta_{b}(2m_{3}+A_{1}^{\prime}-\gamma_{a}\theta_{aa})$
and
$D=\theta_{aa}\bigl{\\{}(B_{1}^{\prime}\theta_{a}+m_{2})A_{1}^{\prime}+\theta_{a}\gamma_{b}\theta_{b}(m_{3}+A_{1}^{\prime})\bigr{\\}}.$
The nature of the curve $f$ can be observe from (16) as

  1. (i)

$f(0)>0$ because $D>0,$

  2. (ii)

$f\to-\infty\,(\infty)$ as $p_{a}\to\infty\,(-\infty)$ because $A<0,$

  3. (iii)

The number of positive real roots of equation (16) can be seen using
Descarte’s rule of sign, mentioned in Table (3).

Therefore, $f(p_{a})=0$ will have either unique or three positive real roots.
Now, to discard the existence of three positive real roots, we use the
nullcline method. From the system (14), we have

$\displaystyle
F_{1}(p_{a},p_{b}):=\frac{m_{a}k_{a}\theta_{b}p_{a}}{\delta_{a}(\theta_{aa}+p_{a})(\theta_{b}+p_{b})}+A_{1}^{\prime}-\gamma_{a}p_{a}$
(17) $\displaystyle
F_{2}(p_{a},p_{b}):=\frac{m_{b}k_{b}\theta_{a}}{\delta_{b}(\theta_{a}+p_{a})}+B_{1}^{\prime}-\gamma_{b}p_{b},$

First to observe the nature of curve $F_{2}(p_{a},p_{b})=0,$ we reduce it in
the form of

$p_{b}=g_{1}(p_{a}):=\frac{1}{k_{2}}\Big{\\{}\frac{m_{b}\theta_{a}}{(\theta_{a}+p_{a})}+B_{1}\Big{\\}}.$
($\mathscr{A}_{2}$)

The following conclusion can be made from ($\mathscr{A}_{2}$)

  1. $\mathscr{A}_{2}$(i)

$g_{1}(0)=\frac{1}{k_{2}}(m_{b}+B_{1})>\frac{B_{1}}{k_{2}}>0$

  2. $\mathscr{A}_{2}$(ii)

$g_{1}(p_{a})\to\frac{B_{1}}{k_{2}}$ as $p_{a}\to\pm\infty$

  3. $\mathscr{A}_{2}$(iii)

$g^{\prime}_{1}(p_{a})<0$, i.e, $g_{1}(p_{a})$ is a decreasing function of
$p_{a}.$

The graph of $F_{2}(p_{a},p_{b})=0$ (or $p_{b}=g_{1}(p_{a})$) is shown in Fig.
9.
Now, to observe the nature of curve $F_{1}(p_{a},p_{b})=0,$ we reduce it in
the form of

$p_{b}=g_{2}(p_{a}):=\frac{-\theta_{b}\left(m_{a}\left(2p_{a}+\theta_{aa}\right)+\left(A_{1}^{\prime}-k_{1}p_{a}\right)\left(\theta_{aa}+p_{a}\right)\right)}{m_{1}p_{a}+\left(A_{1}^{\prime}-k_{1}p_{a}\right)\left(\theta_{aa}+p_{a}\right)}$
($\mathscr{D}$)

The following conclusion can be made from ($\mathscr{D}$)

  1. $\mathscr{D}$(i)

The $p_{a}-$intercept of $g_{2}(p_{a})$ are
$x_{1}=\frac{(2m_{a}+A_{1}-k_{1}\theta_{aa})+\sqrt{(2m_{a}+A_{1}-k_{1}\theta_{aa})^{2}+4k_{1}A_{1}\theta_{aa}}}{2k_{1}}$
and
$x_{2}=\frac{(2m_{a}+A_{1}-k_{1}\theta_{aa})-\sqrt{(2m_{a}+A_{1}-k_{1}\theta_{aa})^{2}+4k_{1}A_{1}\theta_{aa}}}{2k_{1}}$.
We note that $x_{1}>0$ and $x_{2}<0,$ which means we have one positive and one
negative $p_{a}-$intercept. Also $p_{b}-$intercept is negative, i.e.
$g_{2}(0)<0.$

  2. $\mathscr{D}$(ii)

$g_{2}(p_{a})$ is a rational function, and it has two vertical asymptotes at
$x_{3}=\frac{(m_{a}+A_{1}-k_{1}\theta_{aa})+\sqrt{(m_{a}+A_{1}-k_{1}\theta_{aa})^{+}4k_{1}A_{1}\theta_{aa}}}{2k_{1}}$
and
$x_{4}=\frac{(m_{a}+A_{1}-k_{1}\theta_{aa})-\sqrt{(m_{a}+A_{1}-k_{1}\theta_{aa})^{2}+4k_{1}A_{1}\theta_{aa}}}{2k_{1}}$.
We note that $x_{3}>0$ and $x_{4}<0,$ which means we have one positive and one
negative vertical asymptotes. Moreover, $p_{b}=-\theta_{b}$ is a horizontal
asymptote.

We note that $x_{3}<x_{1}.$ For $p_{a}\in(0,x_{3}),$ we have $g_{2}(p_{a})<0$
and for $p_{a}\in(x_{3},x_{1}),$ we have $g_{2}(p_{a})>0.$ Also,
$g_{2}^{\prime}(p_{a})<0$ which implies that $g_{2}(p_{a})$ is strictly
decreasing. The corresponding diagram is shown in Fig. 13 (a). Therefore, from
the above discussion regarding the nature of $g_{1}$ and $g_{2}$, we conclude
that there exists a unique steady state, shown in Fig. 13 (b), and we discard
the case of three positive steady states.

Figure 13: For $n_{a}=n_{b}=n_{aa}=1;$ (a) The graph depicting unique positive
$p_{a}$-intercept of $g_{2}(pa)$ (b) Existence of unique positive steady state
corresponding to Fig. 13 (a).

Case 2: When $n_{a}=n_{b}=1,\,n_{aa}=2,$ then from (15), we get a quartic
polynomial in $p_{a}$

$f(p_{a})=Ap_{a}^{4}+Bp_{a}^{3}+Cp_{a}^{2}+Dp_{a}+E,$ (18)

where
$A=-\gamma_{a}(B_{1}^{\prime}+\gamma_{b}\theta_{b}),\,B=-\gamma_{a}(B_{1}^{\prime}\theta_{a}+m_{2})+B_{1}^{\prime}(m_{3}+A_{1}^{\prime})-\gamma_{a}\gamma_{b}\theta_{a}\theta_{b}+\theta_{b}\gamma_{b}(2m_{3}+A_{1}^{\prime}),\,C=(B_{1}^{\prime}\theta_{a}+m_{2})(m_{3}+A_{1}^{\prime})+\theta_{b}\gamma_{b}\theta_{a}(2m_{3}+A_{1}^{\prime})-\gamma_{a}\theta_{aa}^{2}(B_{1}^{\prime}+\gamma_{b}\theta_{b}),\,D=\theta_{aa}\bigl{\\{}B_{1}^{\prime}A_{1}^{\prime}+\gamma_{b}\theta_{b}(m_{3}+A_{1}^{\prime})\bigr{\\}}-\gamma_{a}\theta_{aa}^{2}(B_{1}^{\prime}\theta_{a}+m_{2}+\gamma_{b}\theta_{b}\theta_{a})$
and
$E={\theta_{aa}}^{2}\bigl{\\{}(B_{1}^{\prime}\theta_{a}+m_{2})A_{1}^{\prime}+\gamma_{b}\theta_{b}\theta_{a}(m_{3}+A_{1}^{\prime})\bigr{\\}}.$
The number of positive real roots of equation (18) can be seen using
Descarte’s rule of sign (mentioned in Table (2)); therefore, $f(p_{a})=0$ will
have either unique or three positive real roots.
Case 3: When $n_{a}=2,\,n_{b}=3,\,n_{aa}=2,$ then from (15), we get a ninth
degree polynomial in $p_{a}$

$f(p_{a})=Ap_{a}^{9}+Bp_{a}^{8}+Cp_{a}^{7}+Dp_{a}^{6}+Ep_{a}^{5}+Fp_{a}^{4}+Gp_{a}^{3}+Hp_{a}^{2}+Ip_{a}+J,$
(19)

where
$A=-\gamma_{a}\,\left({B_{1}^{\prime}}^{3}+{\gamma_{b}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\right),\,B={B_{1}^{\prime}}^{3}\,\left(A_{1}^{\prime}+m_{3}\right)+{\gamma_{b}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\,\left(A_{1}^{\prime}+2\,m_{3}\right),$
$C=-\gamma_{a}\,\left({B_{1}^{\prime}}^{2}\,\left(3\,m_{2}+3\,B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}+B_{1}^{\prime}\,{\theta_{aa}}^{n_{aa}}\right)+{\gamma_{b}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\,\left(3\,{\theta_{a}}^{n_{a}}+{\theta_{aa}}^{n_{aa}}\right)\right),$
$D={B_{1}^{\prime}}^{2}\,\left(\left(A_{1}^{\prime}+m_{3}\right)\,\left(3\,m_{2}+3\,B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)+A_{1}^{\prime}\,B_{1}^{\prime}\,{\theta_{aa}}^{n_{aa}}\right)+{\gamma_{b}}^{n_{b}}\,{\theta_{b}}^{n_{b}}\\\
\left(3\,{\theta_{a}}^{n_{a}}\,\left(A_{1}^{\prime}+2\,m_{3}\right)+{\theta_{aa}}^{n_{aa}}\,\left(A_{1}^{\prime}+m_{3}\right)\right),$
$E=-3\,B_{1}^{\prime}\,\gamma_{a}\,\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)\,\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}+B_{1}^{\prime}\,{\theta_{aa}}^{n_{aa}}\right)-3\,\gamma_{a}\,{\gamma_{b}}^{n_{b}}\,{\theta_{a}}^{n_{a}}\,{\theta_{b}}^{n_{b}}\,\left({\theta_{a}}^{n_{a}}+{\theta_{aa}}^{n_{aa}}\right),$
$F=3\,B_{1}^{\prime}\,\left(\left(A_{1}^{\prime}+m_{3}\right)\,\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)+A_{1}^{\prime}\,B_{1}^{\prime}\,{\theta_{aa}}^{n_{aa}}\right)\,\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)+3\,{\gamma_{b}}^{n_{b}}\,{\theta_{a}}^{n_{a}}\,{\theta_{b}}^{n_{b}}\\\
\left({\theta_{a}}^{n_{a}}\,\left(A_{1}^{\prime}+2\,m_{3}\right)+{\theta_{aa}}^{n_{aa}}\,\left(A_{1}^{\prime}+m_{3}\right)\right),$
$G=-\gamma_{a}\,{\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)}^{2}\,\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}+3\,B_{1}^{\prime}\,{\theta_{aa}}^{n_{aa}}\right)-\gamma_{a}\,{\gamma_{b}}^{n_{b}}\,{\theta_{a}}^{2\,n_{a}}\,{\theta_{b}}^{n_{b}}\,\left({\theta_{a}}^{n_{a}}+3\,{\theta_{aa}}^{n_{aa}}\right),$
$H=\left(\left(A_{1}^{\prime}+m_{3}\right)\,\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)+3\,A_{1}^{\prime}\,B_{1}^{\prime}\,{\theta_{aa}}^{n_{aa}}\right)\,{\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)}^{2}+{\gamma_{b}}^{n_{b}}\,{\theta_{a}}^{2\,n_{a}}\,{\theta_{b}}^{n_{b}}\\\
\left({\theta_{a}}^{n_{a}}\,\left(A_{1}^{\prime}+2\,m_{3}\right)+3\,{\theta_{aa}}^{n_{aa}}\,\left(A_{1}^{\prime}+m_{3}\right)\right),$
$I=-\gamma_{a}\,{\theta_{aa}}^{n_{aa}}\,\left({\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)}^{3}+{\gamma_{b}}^{n_{b}}\,{\theta_{a}}^{3\,n_{a}}\,{\theta_{b}}^{n_{b}}\right),$
and
$J={\theta_{aa}}^{n_{aa}}\,\left(A_{1}^{\prime}\,{\left(m_{2}+B_{1}^{\prime}\,{\theta_{a}}^{n_{a}}\right)}^{3}+{\gamma_{b}}^{n_{b}}\,{\theta_{a}}^{3\,n_{a}}\,{\theta_{b}}^{n_{b}}\,\left(A_{1}^{\prime}+m_{3}\right)\right).$
The number of positive real roots of equation (19) can be obtained using
Descarte’s rule of sign. Due to complexity of polynomial, we discuss the
existence of steady states via numerical simulation.

####  5.1.1 Local stability of steady state

######  Theorem 3.

The steady state $S_{1}^{*}$ is locally asymptotically stable if
$\epsilon_{1}>0,\,\epsilon_{3}>0,\,\epsilon_{4}>0,$ and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}>0,$
where $\epsilon_{i},\,i=1,..,4$ are provided in the proof.

###### Proof.

At $S_{1}^{*},$ the Jacobian matrix of the model system (14) is given by

$J_{S_{1}^{*}}=\left[\begin{array}[]{cccc}-\gamma_{a}&0&X&Y\\\
0&-\gamma_{b}&Z&0\\\ k_{a}&0&-\delta_{a}&0\\\
0&-k_{b}&0&-\delta_{b}\end{array}\right].$

where
$X=\frac{\mathrm{m}_{a}\,n_{aa}\,{\theta_{aa}}^{n_{aa}}\,p_{a}^{n_{aa-1}}}{{\left({\theta_{aa}}^{n_{aa}}+p_{a}^{n_{aa}}\right)}^{2}},\,Y=-\frac{\mathrm{m}_{a}\,n_{b}\,{\theta_{b}}^{n_{b}}\,p_{b}^{n_{b}-1}}{{\left({\theta_{b}}^{n_{b}}+p_{b}^{n_{b}}\right)}^{2}},$
and
$Z=-\frac{\mathrm{m}_{b}\,n_{a}\,{\mathrm{p_{a}}}^{n_{a}-1}\,{\theta_{a}}^{n_{a}}}{{\left({\mathrm{p_{a}}}^{n_{a}}+{\theta_{a}}^{n_{a}}\right)}^{2}}.$
The characteristic equation of $J_{S_{1}^{*}}$ is given as

$\lambda^{4}+\epsilon_{1}\lambda^{3}+\epsilon_{2}\lambda^{2}+\epsilon_{3}\lambda+\epsilon_{4}=0.$
(20)

Here,
$\epsilon_{1}=\gamma_{a}+\gamma_{b}+\delta_{a}+\delta_{b},\epsilon_{2}=\gamma_{a}\gamma_{b}+(\gamma_{a}+\gamma_{b})(\delta_{a}+\delta_{b})+\delta_{a}\delta_{b}-k_{a}X,\,\epsilon_{3}=\delta_{a}\delta_{b}(\gamma_{a}+\gamma_{b})+(\delta_{a}+\delta_{b})(\gamma_{a}\gamma_{b})-k_{a}(\delta_{b}+\gamma_{b})X,$
and $\epsilon_{4}=\gamma_{a}\gamma_{b}\delta_{a}\delta_{b}-k_{a}k_{b}YZ-
k_{a}\gamma_{b}\delta_{b}X.$ By the Routh-Hurwitz criteria, all the roots of
(20) have negative real parts provided
$\epsilon_{1}>0,\,\epsilon_{3}>0,\,\epsilon_{4}>0,$ and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}>0.$
∎

####  5.1.2 Saddle Node Bifurcation

We find the transversality condition the same as discussed in Subsection
3.1.2. Define $f=(f_{1},f_{2},f_{3},f_{4})^{T},$ where $f_{1}$ are defined as

$f_{1}=m_{a}\biggl{(}\frac{(\frac{p_{b}}{\theta_{b}})^{n_{b}}}{1+(\frac{p_{b}}{\theta_{b}})^{n_{b}}}+\frac{1}{1+(\frac{p_{b}}{\theta_{b}})^{n_{b}}}\biggr{)}-\gamma_{a}r_{a}+A_{1},$

and the expression of $f_{2},\,f_{3},$ and $f_{4}$ are same as Subsection
3.1.2. Then the model system (14) experiences a saddle-node bifurcation at
$(S_{SN}^{*},\theta_{b}^{SN})$, if the following transversality conditions
[36] are satisfied:

$w^{T}f_{\theta_{b}}(S_{SN}^{*},\theta_{b}^{SN})=w_{1}\frac{m_{a}n_{b}{\theta_{b}^{n_{b}-1}}^{SN}{p_{b}^{n_{b}}}^{{SN}^{*}}}{(\mathcal{X}_{2})^{2}}\neq
0.$ ($\mathcal{T}_{33}$)

and

$\displaystyle w^{T}[D^{2}f(S_{SN}^{*},\theta_{b}^{SN})(v,v)]=$ $\displaystyle
w_{1}\Biggl{\\{}\frac{\mathcal{O}v_{3}^{2}}{\mathcal{X}_{1}^{3}}\bigg{(}\theta_{aa}^{n_{aa}}(n_{aa}-1){p_{a}^{(n_{aa}-2)}}^{{SN}^{*}}-(n_{aa}+1){p_{a}^{(2n_{aa}-2)}}^{{SN}^{*}}\bigg{)}$
$\displaystyle-\frac{v_{4}^{2}\mathcal{O}_{1}}{\mathcal{X}_{2}^{3}}\bigg{(}{\theta_{b}^{n_{b}}}^{SN}(n_{b}-1){p_{b}^{(n_{b}-2)}}^{{SN}^{*}}-(n_{b}+1){p_{b}^{(2n_{b}-2)}}^{{SN}^{*}}\bigg{)}\Biggr{\\}}$
$\displaystyle-\frac{w_{2}v_{3}^{2}m_{b}n_{a}\theta_{a}^{n_{a}}}{(\theta_{a}^{n_{a}}+{p_{a}^{n_{a}}}^{{SN}^{*}})^{3}}\bigg{\\{}{p_{a}^{n_{a}-2}}^{{SN}^{*}}\theta_{a}^{n_{a}}(n_{a}-1)-(n_{a}+1){p_{a}^{2n_{a}-2}}^{{SN}^{*}}\bigg{\\}}$
$\displaystyle\neq 0,$ ($\mathcal{T}_{44}$)

where
$\mathcal{X}_{1}=\theta_{aa}^{n_{aa}}+{p_{a}^{n_{aa}}}^{{SN}^{*}},\,\mathcal{X}_{2}=\theta_{b}^{n_{b}}+{p_{b}^{n_{b}}}^{{SN}^{*}},\,\mathcal{O}=n_{aa}m_{a}\theta_{aa}^{n_{aa}},$
and $\mathcal{O}_{1}=m_{a}n_{b}\theta_{b}^{n_{b}}.$

####  5.1.3 Numerical simulation and discussion corresponding to model system
(14)

The purpose of this section is to provide validation for the theoretical
results that were covered in Subsection 5.1. We have talked about the system’s
behavior in terms of monostability, bistability, and tristability with
fictitious parameters.

######  Example 7.

$n_{a}=n_{b}=n_{aa}=1$ (case 1) and the other parameters are chosen as
$\delta_{a}=0.09,\,\delta_{b}=0.9,\,\gamma_{a}=1,\,\gamma_{b}=0.9,\,k_{a}=0.25,\,k_{b}=0.9,\,m_{b}=0.1,\,m_{a}=2,\,\theta_{a}=0.4,\,A_{1}=0.0035,\,B_{1}=7.85,\,\theta_{aa}=20,\,,\,\theta_{b}=(0.00001,50).$
For this set of parameters, we obtain a unique stable steady state as
discussed in Case 1. The corresponding bifurcation plot is shown in Fig. 14
(a). To ensure the stability of $S_{1}^{*},$ we choose $\theta_{b}=10$ from
the range of $\theta_{b}.$ The coefficients of characteristic equation at
$S_{1}^{*}$ are calculated as
$\epsilon_{1}=2.89>0,\,\epsilon_{2}=2.8445>0,\,\epsilon_{3}=1.0133>0,\,\epsilon_{4}=0.0790>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.0587>0.$
Therefore, from Theorem 3, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
$\theta_{b}$, shown in Fig. 14 (b). We can see that all trajectories starting
from any initial points are converging to
$S_{1}^{*}=(1.39595,8.73261,3.87764,8.73261).$

Figure 14: For $n_{a}=1,\,n_{b}=1,\,n_{aa}=1$; (a) Bifurcation plot with the
unique stable steady state (b) Phase portrait at $\theta_{b}=10.$

######  Example 8.

$n_{a}=n_{b}=1,\,n_{aa}=2$ (case 2), and the other parameters are chosen as:
$\delta_{a}=0.363,\,\delta_{b}=0.35,\,\gamma_{a}=0.84,\,\gamma_{b}=0.9,\,k_{a}=2,\,k_{b}=6.3,\,m_{b}=5.7,\,m_{a}=3.4,\,\theta_{a}=6.2,\,A_{1}=0.1,\,B_{1}=0.01,\,\theta_{aa}=13.1,\,\theta_{b}=(1.6,6.5).$
In this range of $\theta_{b},$ there is unique ($S_{1}^{*}$) steady state in
$\theta_{b}^{1}$, three ($S_{1}^{*},\,S_{2}^{*},\,S_{3}^{*}$) in
$\theta_{b}^{2}$ and unique ($S_{3}^{*}$) in $\theta_{b}^{3},$ such that
$\theta_{b}=\theta_{b}^{1}\cup\theta_{b}^{2}\cup\theta_{b}^{2}$ where
$\theta_{b}^{1}=(1.6,{\theta_{b}^{SN}}^{1}),\,\theta_{b}^{2}=[{\theta_{b}^{SN}}^{1},{\theta_{b}^{SN}}^{2}]$
and $\theta_{b}^{2}=({\theta_{b}^{SN}}^{2},6.5).$ The steady states
$S_{1}^{*},\,S_{3}^{*}$ are stable and ($S_{2}^{*}$) is unstable. The
correspond bifurcation diagram is shown in Fig. 15 (a).
To ensure the stability of $S_{1}^{*},$ we chosen
$\theta_{b}=1.89\in\theta_{b}^{1}$ and the coefficients of characteristic
equation at $S_{1}^{*}=(0.239206,5.23418,1.31794,94.2152)$ are calculated as
$\epsilon_{1}=2.453>0,\,\epsilon_{2}=2.0213>0,\,\epsilon_{3}=0.6321>0,\,\epsilon_{4}=0.0583>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=2.3838>0.$
Therefore, from Theorem 3, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
parameters, shown in Fig. 15 (b). We can see that all trajectories starting
from any initial points are converging to $S_{1}^{*}.$
Further, we chosen $\theta_{b}=3.81\in\theta_{b}^{2}$ and the coefficients of
characteristic equation at $S_{1}^{*}=(0.416123,4.6347,2.29269,83.4245)$ are
calculated as
$\epsilon_{1}=2.453>0,\,\epsilon_{2}=1.9526>0,\,\epsilon_{3}=0.5463>0,\,\epsilon_{4}=0.0317>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=2.1276>0.$
The coefficients of characteristic equation at
$S_{2}^{*}=(1.1931,3.08517,6.57356,55.5331)$ are calculated as
$\epsilon_{1}=2.453>0,\,\epsilon_{2}=1.7912>0,\,\epsilon_{3}=0.3445<0,\,\epsilon_{4}=-0.0187<0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=1.5077>0.$
The coefficients of characteristic equation at
$S_{3}^{*}=(3.15301,1.67693,17.3719,30.1848)$ are calculated as
$\epsilon_{1}=2.453>0,\,\epsilon_{2}=1.9428>0,\,\epsilon_{3}=0.5339>0,\,\epsilon_{4}=0.03001>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=2.0785>0.$
Therefore, from Theorem 3, we conclude that $S_{1}^{*},\,S_{3}^{*}$ is locally
asymptotically stable and $S_{2}^{*}$ is unstable. We draw the phase portrait
for the considered parameters, shown in Fig. 15 (c). We can see that all
trajectories are converging to $S_{1}^{*},\,S_{3}^{*},$ and nearby
trajectories of $S_{2}^{*}$ are moving away from $S_{2}^{*}.$ This shows the
bistable behavior of the model system (14).
We chosen $\theta_{b}=5.94\in\theta_{b}^{3},$ the coefficients of
characteristic equation at $S_{3}^{*}=(3.7429,\\\ 1.47508,20.622,26.5515)$ are
calculated as
$\epsilon_{1}=2.453>0,\,\epsilon_{2}=1.9886>0,\,\epsilon_{3}=0.5912>0,\,\epsilon_{4}=0.0417>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=2.2838>0.$
Therefore, from Theorem 3, we conclude that $S_{3}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the chosen parameters,
shown in Fig. 15 (d). We can see that all trajectories starting from any
initial points are converging to $S_{3}^{*}.$
In this example, we also discuss the case of saddle-node bifurcation. We have
seen that the system (14) has two steady states $S_{2}^{*},\,S_{3}^{*}$ in
$\theta_{b}^{2}$ such that as the value of $\theta_{b}$ decreases, the two
steady states collide at $\theta_{b}^{{{SN}^{1}}}=2.154867833473525$ and
denoted as $S_{{SN}^{1}}^{*}=(2.12792,2.20183,11.7241,39.6329).$ The Jacobian
matrix $J=Df(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})$ has a simple
eigenvalue $\lambda=0$ and $v=(-0.0748,0.0504,-0.4121,0.9067)^{T},$
$w=(-0.9192,0.0771,-0.3861,0.0110)^{T}$, are the eigenvectors of $J$ and
$J^{T},$ respectively. Both the tranversality conditions
$w^{T}f_{\theta_{b}}(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})=-0.0709\neq 0$
and
$w^{T}[D^{2}f(S_{{SN}^{1}}^{*},\theta_{b}^{{{SN}^{1}}})(v,v)]=0.0015\neq 0$
are satisfied. Hence, the system (14) experiences saddle-node bifurcation at
$\theta_{b}=\theta_{b}^{{{SN}^{1}}}.$ Also, the system (14) has two steady
states $S_{1}^{*},\,S_{2}^{*}$ in $\theta_{b}^{2}$ such that as the value of
$\theta_{b}$ increases, the two steady states collide at
$\theta_{b}^{{{SN}^{2}}}=4.804927166491674$ and denoted as
$S_{{SN}^{2}}^{*}=(0.7238,3.8654,3.9879,69.5764).$ The Jacobian matrix
$J=Df(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})$ has a simple eigenvalue
$\lambda=0$ and
$v=(0.0263,-0.0549,0.1450,-0.9876)^{T},\,w=(0.9206,-0.0544,0.3866,-0.0078)^{T}$,
are the eigenvectors of $J$ and $J^{T},$ respectively. Both the tranversality
conditions
$w^{T}f_{\theta_{b}}(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})=0.0394\neq 0$
and
$w^{T}[D^{2}f(S_{{SN}^{2}}^{*},\theta_{b}^{{{SN}^{2}}})(v,v)]=0.0004194\neq 0$
are satisfied. Hence, the system (14) experiences saddle-node bifurcation at
$\theta_{b}=\theta_{b}^{{{SN}^{2}}}.$
The system (14) exhibits a hysteresis effect in $\theta_{b}^{2}$ where
multiple steady states coexist, as shown in Fig 15 (a). The two outer steady
states are stable, while the interior steady state (red) is unstable.
Again, in this case, the parameters are chosen as
$\delta_{a}=0.09,\,\delta_{b}=0.9,\,\gamma_{a}=1,\,\gamma_{b}=0.9,\,k_{a}=0.25,\,k_{b}=0.9,\,m_{b}=0.1,\,m_{a}=2,\,\theta_{a}=0.4,\,A_{1}=0.0035,\,B_{1}=7.85,\,\theta_{aa}=20,\,\theta_{b}=(0.00001,50).$
For this set of parameters, we obtain a unique stable steady state as
discussed in Case 2 (OR2). The corresponding bifurcation plot is shown in Fig.
15 (e). To ensure the stability of $S_{1}^{*},$ we choose $\theta_{b}=10$ from
the range of $\theta_{b}.$ The coefficients of characteristic equation at
$S_{1}^{*}$ are calculated as
$\epsilon_{1}=2.89>0,\,\epsilon_{2}=2.8546>0,\,\epsilon_{3}=1.0316>0,\,\epsilon_{4}=0.0669>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=6.88757>0.$
Therefore, from Theorem 3, we conclude that $S_{1}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the considered
$\theta_{b}$, shown in Fig. 15 (f). We can see that all trajectories starting
from any initial points are converging to
$S_{1}^{*}=(1.11812,8.7349,3.1059,8.7349).$

Figure 15: For $n_{a}=1,\,n_{b}=1,\,n_{aa}=2$; (a) Bifurcation plot exhibiting
hysteresis effect with three steady states (b) Phase portrait at
$\theta_{b}=1.89$ (c) Phase portrait at $\theta_{b}=3.81$ depicting
bistability (d) Phase portrait at $\theta_{b}=5.94$ (e) Bifurcation plot with
the unique stable steady state (f) Phase portrait at $\theta_{b}=10.$

######  Example 9.

$n_{a}=2,\,n_{b}=3,\,n_{aa}=2$ (case 3) and the other parameters are chosen as
$\delta_{a}=0.32,\,\delta_{b}=0.35,\,\gamma_{a}=0.25,\,\gamma_{b}=1,\,k_{a}=0.5,\,k_{b}=0.91,\,m_{b}=15.5,\,m_{a}=3.4,\,\theta_{a}=9.4,\,A_{1}=0.2,\,B_{1}=0.01,\,\theta_{aa}=11.67,\,\theta_{b}\in(3,11).$
In this range of $\theta_{b},$ there is three
($S_{1}^{*},\,S_{2}^{*},\,S_{3}^{*}$) steady states in $\theta_{b}^{1}$, five
($S_{1}^{*},\,S_{2}^{*},\,S_{3}^{*},\,S_{4}^{*},\,S_{5}*$) in
$\theta_{b}^{2},$ three ($S_{1}^{*},\,S_{2}^{*},\,S_{5}*)$ in
$\theta_{b}^{3},$ and unique ($S_{5}^{*}$) in $\theta_{b}^{4},$ such that
$\theta_{b}=\theta_{b}^{1}\cup\theta_{b}^{2}\cup\theta_{b}^{2}\cup\theta_{b}^{4}$
where
$\theta_{b}^{1}=(3,3.51),\,\theta_{b}^{2}=(3.52,3.95),\,,\,\theta_{b}^{3}=(3.96,10.21)$
and $\theta_{b}^{4}=(10.22,11).$ The steady states
$S_{1}^{*},\,S_{3}^{*},\,S_{5}^{*}$ are stable and $S_{2}^{*},\,S_{4}^{*}$ are
unstable. The corresponding bifurcation diagram is shown in Fig. 16 (a).
To ensure the tristability in $\theta_{b}^{2},$ we chosen
$\theta_{b}=3.7\in\theta_{b}^{2}$ and the coefficients of characteristic
equation at $S_{1}^{*}=(1.10066,15.008,1.71979,39.0207)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.1584>0,\,\epsilon_{3}=0.2520>0,\,\epsilon_{4}=0.0135>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.4470>0.$
The coefficients of characteristic equation at
$S_{2}^{*}=(4.69214,9.64743,7.33148,25.0833)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.1054>0,\,\epsilon_{3}=0.1805>0,\,\epsilon_{4}=-0.0055<0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.3708>0.$
The coefficients of characteristic equation at
$S_{3}^{*}=(10.4286,3.88023,16.2947,10.0886)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.1527>0,\,\epsilon_{3}=0.2444>0,\,\epsilon_{4}=0.0043>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.4653>0.$
The coefficients of characteristic equation at
$S_{4}^{*}=(16.8118,1.7695,26.2684,4.60071)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.1817>0,\,\epsilon_{3}=0.2834>0,\,\epsilon_{4}=-0.0052<0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.5819>0.$
The coefficients of characteristic equation at
$S_{5}^{*}=(23.5993,0.95581,36.874,2.48511)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.1919>0,\,\epsilon_{3}=0.2972>0,\,\epsilon_{4}=0.0093>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.5576>0.$
Therefore, from Theorem 3, we conclude that
$S_{1}^{*},\,S_{3}^{*},\,S_{5}^{*}$ are locally asymptotically stable and
$S_{2}^{*},\,S_{4}^{*}$ are unstable. We draw the phase portrait for the
considered parameters, shown in Fig. 16 (b). We can see that all trajectories
starting from any initial points are converging to
$S_{1}^{*},\,S_{3}^{*},\,S_{5}^{*}$ and nearby trajectories of
$S_{2}^{*},\,S_{4}^{*}$ are moving away from $S_{2}^{*},\,S_{4}^{*}.$ This
shows the tristable behavior of the model system (14).
Further, to ensure the bistability in $\theta_{b}^{3},$ we chosen
$\theta_{b}=7\in\theta_{b}^{3},$ and the coefficients of characteristic
equation at $S_{1}^{*}=(1.25204,14.8665,1.95631,38.653)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.1533>0,\,\epsilon_{3}=0.2451>0,\,\epsilon_{4}=0.0114>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.4407>0.$
The coefficients of characteristic equation at
$S_{2}^{*}=(3.84503,11.0147,6.00785,\\\ 28.6382)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.1058>0,\,\epsilon_{3}=0.1810<0,\,\epsilon_{4}=-0.0072<0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.3782>0.$
The coefficients of the characteristic equation at
$S_{5}^{*}=(26.7174,0.757962,41.7459,1.9707)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.194>0,\,\epsilon_{3}=0.3001>0,\,\epsilon_{4}=0.0244>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.5081>0.$
Therefore, from Theorem 3, we conclude that $S_{1}^{*},\,S_{5}^{*}$ are
locally asymptotically stable and $S_{2}^{*}$ is unstable. We draw the phase
portrait for the considered parameters, shown in Fig. 16 (c). We can see that
all trajectories are converging to $S_{1}^{*},\,S_{5}^{*},$ and nearby
trajectories of $S_{2}^{*}$ are moving away from $S_{2}^{*}.$ This shows the
bistable behavior of the model system (14).
To ensure the mono-stability in $\theta_{b}^{4},$ we chosen
$\theta_{b}=10.5\in\theta_{b}^{4},$ and the coefficients of characteristic
equation at $S_{5}^{*}=(26.9444,0.746008,42.1006,1.93962)$ are calculated as
$\epsilon_{1}=1.92>0,\,\epsilon_{2}=1.1941>0,\,\epsilon_{3}=0.3003>0,\,\epsilon_{4}=0.0256>0,$
and
$\epsilon_{1}\epsilon_{2}\epsilon_{3}-\epsilon_{1}^{2}\epsilon_{4}-\epsilon_{3}^{2}=0.5038>0.$
Therefore, from Theorem 3, we conclude that $S_{5}^{*}$ is locally
asymptotically stable. We draw the phase portrait for the chosen parameters,
shown in Fig. 16 (d). We can see that all trajectories starting from any
initial points are converging to $S_{5}^{*}.$

Figure 16: For $n_{a}=2,\,n_{b}=3,\,n_{aa}=2$; (a) Bifurcation plot with five
steady states (b) Phase portrait at $\theta_{b}=3.7$ depicting tristability
(c) Phase portrait at $\theta_{b}=7$ depicting bistability (d) Phase portrait
at $\theta_{b}=10.5.$

###  5.2 Conclusion and biological interpretation for this section

In this section, we consider a non-competitive type-II $\mathcal{OR}$ logic as
an alternate to competitive OR logic. Here we assume that the two genes $a$
and $b$ don’t compete for binding and any binding events can elicit gene
expression. Our results overlap with that of AND logic by showing that
bistability requires dimeric (Hill coefficients of degree 2) autoloop and cell
state transition can occur only through abrupt transitions (hysteresis).
Nevertheless, the main result of this section is tristability, which we show
is possible only in higher degrees of multimerization (e.g., ($n_{a},\
n_{aa},\ n_{b})\ =\ (2,\ 2,\ 3$)). One manifestation of tristability is
phenotypic heterogeneity (concurrent occurrence of three phenotypes) and we
also observe that cell transition occurs abruptly after a certain threshold of
$\theta_{b}$ is crossed (Fig. 16a). We thus conclude that an autoregulated
two-node positive feedback circuit with non-competitive OR logic (type-II
$\mathcal{OR}$) is the ”minimal” circuit to exhibit multistability under
higher degrees of multimerization.

##  6 Impact of basal production rate on model dynamics

Basal production rate or leakage defines the production of mRNA in the absence
of any activatory or inhibitory interacting transcription factors. For
instance, the basal production of $r_{a}$ will be the rate of production of
$r_{a}$ in ”normal” conditions, i.e., when the self-activation due to $a$
itself or inhibition by $b$ is absent. While modeling gene expression
(transcription) processes, leakage is sometimes neglected at the cost of
simplifying modeling efforts and analytical calculations. However, in all of
our model systems, (1), (8) and (14), we have considered the leakage rates of
mRNA’s, $r_{a}$ and $r_{b}$ as $A_{1}$ and $B_{1}$, respectively.

We wanted to investigate whether the absence of the leakage rate of $r_{a}$
will have any impact on the model dynamics. To understand that, we compared
the mRNA concentration $r_{a}$ for the model systems (1), (8) and (14) for the
zero and non-zero values of $A_{1},$ as shown in Fig. 17 (a) and (b),
respectively. From Fig. 17, we observed that the level of mRNA concentration
$r_{a}$ is high for the model with type-II $\mathcal{OR}$ logic (model system
14) (blue color curve) and low for the model with type-I OR logic (model
system 1) (orange color curve). Also, it is evident that the level of mRNA
concentration $r_{a}$ is high for all the three model systems with non-zero
$A_{1}$ (Fig. 17b) compared to $A_{1}=0$ case (Fig. 17a). However, the
structure of the time-dynamics curves of mRNA concentration remains largely
unaltered. Our results thus show that leakage only plays the role of ”scaling”
in the model and its incorporation in the model systems won’t have any impact
on the overall qualitative dynamics of the network. We thus conclude that our
results hold true even when the leakage is completely ignored.

Figure 17: Comparative study of temporal dynamics of mRNA concentration
$r_{a}$ of gene $a$ using three regulatory logics: Type-I OR logic, AND logic,
and Type-II $\mathcal{OR}$ logic.(a) when $A_{1}=0$ and (b) when $A_{1}\neq
0$. The parameter values for these plots are same as that for Fig. 4 (f)

.

##  7 Discussion and Final Conclusion

Cell fate switching is a dynamic phenomenon often tied to regulatory network
motifs that at the cellular level define the computational machinery of life.
Most of these network motifs define molecular switches exhibiting diverse
qualitative behaviors such as bistability, catastrophes, and hysteresis [1].
The most prominent examples of molecular switches involving minimal circuitry
are the two-component positive feedback network motifs resulting from mutual
repression or mutual activation of two genes. Traditional studies have shown
that these motifs can exhibit at most bistable dynamics (two stable states or
attractors) allowing the system to alternately switch between two states.
Further, the underlying cause of bistable dynamics is attributed to multimeric
regulation (higher order Hill coefficient), in contrast to monomeric
regulation that gives rise to only monostable dynamics. The biological
programs underlying cell fate decision-making are however not just restricted
to mono-and bistable regimes of dynamics. In fact, higher-order dynamics such
as tristability (three stable states/attractors/phenotypes) is now also
prominently observed in biological mechanisms. This happens in differentiation
programs such as differentiation of naive CD4+ T cells [37] as well as in
diseases-inducing processes such as epithelial-mesenchymal transition (EMT) in
which non-motile epithelial cells switch to mesenchymal and hybrid epithelial-
mesenchymal cell fates which have migratory and invasive traits that often
causes metastasis of carcinomas [26]. Nevertheless, the minimal network motif
that can exhibit tristability, or for that matter, mon, bi-and tri stability,
remains largely unexplored.

In our previous study on networks underlying and driving EMT, we used a
“numerical” approach to report that a positive feedback loop formed by mutual
repression of two genes, miR200 ( microRNA) and ZEB (mRNA), with one gene,
ZEB, self-activated can exhibit tristability. However, this being a microRNA-
mRNA circuit involves translational repression and active degradation of ZEB
due to micro-RNA, making it a complicated mechanism. We asked whether a
generic positive feedback loop formed by mutual repression of two
transcription factors (rather than microRNAs and mRNA’s) with an autoloop, can
be the minimal circuit to exhibit lower as well as higher-order stability
regimes such as mon, bi-and tri stability. While small feedback circuits with
and without delays have been studied previously [38, 39, 40, 41, 42, 43, 44],
we still lack answers with rigorous analytical foundations to the following
crucial questions: (a) How is multimerization associated with the order of
stability? (b) Does the “logic” of regulation have any functional role in
multistability? In this work, we address these questions using a rigorous
analytical (and numerical) approach that allows us to reach the following
conclusions summarised in Table 4.

We show that a positive feedback loop with an autoloop can exhibit diverse
dynamical behaviors that depend on the degree of multimerization as well as
the “logic” of regulation. The autoloop can be integrated into the core
circuit through three types of logic gates: OR (Type-I and Type-II), and AND.
With Type-I OR logic, the monomeric circuit with monomeric autoloop can show
biphasic kinetics as well as hysteresis (bistability). However, with AND and
Type-II logics the circuit can show hyperbolic saturation kinetics - a
manifestation of monostability. We show that the necessary geometrical
constraint for hysteresis (bistability) with AND and Type-II logics imposes
the existence of (at least) dimeric autoloop. We also analyzed higher degrees
of multimerization and reported the prevalence of biphasic as well as mono-and
bistable dynamics. Finally, and strikingly, we show the possibility of
tristability in this circuit. The two requirements for tistable dynamics are:
(a) the degrees of multimerization for the two regulations in the core circuit
should be (at least) two and three, (b) the degree of multimerization of
autoloop should be (at least) two (i. e., dimeric). Taken together, our
analysis shows how different multimerization and logical constraints in this
minimal circuit can work together to produce numerous types of dynamical
scenarios including monostability, hysteresis (bistability), biphasic
kinetics, and importantly, tristability. While hysteresis translates to forced
and abrupt cell state transitions, biphasic kinetics explains $smooth\ cell\
state\ swap$ without a step-like switch. The biological network (i.e.,
miR200-ZEB circuit) corresponding to this generic circuit underlies EMT, and
drives cell fate transition and metastasis in carcinomas [45]. Besides, it has
been shown (experimentally) to exhibit hysteric as well non-hysteric dynamics,
with only the hysteric EMT enabling lung metastasis [46], our results can thus
have crucial implications in furthering our understanding of EMT mechanism due
to miR200-ZEB feedback loop and can provide inputs to control cell-fate
switching during metastasis of carcinomas.

$(n_{a},n_{b},n_{aa})$ | Type-I OR logic | AND logic | Type-II $\mathcal{OR}$ logic
---|---|---|---
Monostable |  | Bistable
---
(Hysteresis)
Tristable | Monostable |  | Bistable
---
(Hysteresis)
Tristable | Monostable |  | Bistable
---
(Hysteresis)
Tristable
(1,1,1) | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | - | - | $\checkmark$ | - | -
(1,1,2) | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | -
(2,1,2) | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | -
(2,2,1) | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | -
(2,2,2) | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | -
(2,3,2) | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | $\checkmark$
(3,3,3) | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | - | $\checkmark$ | $\checkmark$ | $\checkmark$

Table 4: Summary of different dynamical scenarios exhibited by the network.

Author Contributions
MR: Conceptualization, supervision, funding acquisition. AS: Analytical
calculations and numerical simulations. Both authors analysed the data,
discussed the results, and wrote the manuscript.

Acknowledgments
This work is supported by the Department of Science and Technology, India
[Grant No. DST/INSPIRE/04/2020/001492] and the Science and Engineering
Research Board, India [Grant No. CRG/2023/006432] to Mubasher Rashid.

## References

  * [1] M Sáez, J Briscoe, and David A Rand.  Dynamical landscapes of cell fate decisions.  Interface focus, 12(4):20220002, 2022.
  * [2] Joseph X Zhou and Sui Huang.  Understanding gene circuits at cell-fate branch points for rational cell reprogramming.  Trends in genetics, 27(2):55–62, 2011.
  * [3] Laura Prochazka, Yaakov Benenson, and Peter W Zandstra.  Synthetic gene circuits and cellular decision-making in human pluripotent stem cells.  Current Opinion in Systems Biology, 5:93–103, 2017.
  * [4] Gábor Balázsi, Alexander Van Oudenaarden, and James J Collins.  Cellular decision making and biological noise: from microbes to mammals.  Cell, 144(6):910–925, 2011.
  * [5] Raúl Guantes and Juan F Poyatos.  Multistable decision switches for flexible control of epigenetic differentiation.  PLoS computational biology, 4(11):e1000235, 2008.
  * [6] Dasong Huang and Ruiqi Wang.  Exploring the mechanisms of cell reprogramming and transdifferentiation via intercellular communication.  Physical Review E, 102(1):012406, 2020.
  * [7] Shubham Tripathi, Herbert Levine, and Mohit Kumar Jolly.  The physics of cellular decision making during epithelial-mesenchymal transition.  Annual Review of Biophysics, 49:1–18, 2020.
  * [8] Shubham Tripathi, Priyanka Chakraborty, Herbert Levine, and Mohit Kumar Jolly.  A mechanism for epithelial-mesenchymal heterogeneity in a population of cancer cells.  PLoS computational biology, 16(2):e1007619, 2020.
  * [9] Dandan Li, Lingyun Xia, Pan Huang, Zidi Wang, Qiwei Guo, Congcong Huang, Weidong Leng, and Shanshan Qin.  Heterogeneity and plasticity of epithelial–mesenchymal transition (emt) in cancer metastasis: Focusing on partial emt and regulatory mechanisms.  Cell proliferation, 56(6):e13423, 2023.
  * [10] Meredith S Brown, Behnaz Abdollahi, Owen M Wilkins, Hanxu Lu, Priyanka Chakraborty, Nevena B Ognjenovic, Kristen E Muller, Mohit Kumar Jolly, Brock C Christensen, Saeed Hassanpour, et al.  Phenotypic heterogeneity driven by plasticity of the intermediate emt state governs disease progression and metastasis in breast cancer.  Science advances, 8(31):eabj8002, 2022.
  * [11] Mubasher Rashid, Kishore Hari, John Thampi, Nived Krishnan Santhosh, and Mohit Kumar Jolly.  Network topology metrics explaining enrichment of hybrid epithelial/mesenchymal phenotypes in metastasis.  PLOS Computational Biology, 18(11):e1010687, 2022.
  * [12] Javier Macía, Stefanie Widder, and Ricard Solé.  Why are cellular switches boolean? general conditions for multistable genetic circuits.  Journal of theoretical biology, 261(1):126–135, 2009.
  * [13] Marcelline Kaufman, Christophe Soule, and René Thomas.  A new necessary condition on interaction graphs for multistationarity.  Journal of theoretical biology, 248(4):675–685, 2007.
  * [14] David Angeli, James E Ferrell Jr, and Eduardo D Sontag.  Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems.  Proceedings of the National Academy of Sciences, 101(7):1822–1827, 2004.
  * [15] German Enciso and Eduardo D Sontag.  Monotone systems under positive feedback: multistability and a reduction theorem.  Systems & control letters, 54(2):159–168, 2005.
  * [16] David Angeli and Eduardo D Sontag.  Multi-stability in monotone input/output systems.  Systems & Control Letters, 51(3-4):185–202, 2004.
  * [17] Christophe Soulé.  Graphic requirements for multistationarity.  ComPlexUs, 1(3):123–133, 2003.
  * [18] Timothy S Gardner, Charles R Cantor, and James J Collins.  Construction of a genetic toggle switch in escherichia coli.  Nature, 403(6767):339–342, 2000.
  * [19] Qing Zhou, Hiram Chipperfield, Douglas A Melton, and Wing Hung Wong.  A gene regulatory network in mouse embryonic stem cells.  Proceedings of the National Academy of Sciences, 104(42):16438–16443, 2007.
  * [20] Ahmadreza Ghaffarizadeh, Nicholas S Flann, and Gregory J Podgorski.  Multistable switches and their role in cellular differentiation networks.  BMC bioinformatics, 15:1–13, 2014.
  * [21] Jinfang Zhu, Hidehiro Yamane, and William E Paul.  Differentiation of effector cd4 t cell populations.  Annual review of immunology, 28:445–489, 2009.
  * [22] Atchuta Srinivas Duddu, Sarthak Sahoo, Souvadra Hati, Siddharth Jhunjhunwala, and Mohit Kumar Jolly.  Multi-stability in cellular differentiation enabled by a network of three mutually repressing master regulators.  Journal of the Royal Society Interface, 17(170):20200631, 2020.
  * [23] Laurane De Mot, Didier Gonze, Sylvain Bessonnard, Claire Chazaud, Albert Goldbeter, and Geneviève Dupont.  Cell fate specification based on tristability in the inner cell mass of mouse blastocysts.  Biophysical journal, 110(3):710–722, 2016.
  * [24] En Li.  Chromatin modification and epigenetic reprogramming in mammalian development.  Nature Reviews Genetics, 3(9):662–673, 2002.
  * [25] Mingyang Lu, Mohit Kumar Jolly, Ryan Gomoto, Bin Huang, Jose Onuchic, and Eshel Ben-Jacob.  Tristability in cancer-associated microrna-tf chimera toggle switch.  The journal of physical chemistry B, 117(42):13164–13174, 2013\.
  * [26] Mingyang Lu, Mohit Kumar Jolly, Herbert Levine, José N Onuchic, and Eshel Ben-Jacob.  Microrna-based regulation of epithelial–hybrid–mesenchymal fate determination.  Proceedings of the National Academy of Sciences, 110(45):18144–18149, 2013.
  * [27] Mubasher Rashid, Brasanna M Devi, and Malay Banerjee.  Combinatorial cooperativity in mir200-zeb feedback network can control epithelial-mesenchymal transition.  Bulletin of Mathematical Biology, 86(5):1–23, 2024.
  * [28] Vijay Chickarmane, Carl Troein, Ulrike A Nuber, Herbert M Sauro, and Carsten Peterson.  Transcriptional dynamics of the embryonic stem cell switch.  PLoS computational biology, 2(9):e123, 2006.
  * [29] Daniel Huang, William J Holtz, and Michel M Maharbiz.  A genetic bistable switch utilizing nonlinear protein degradation.  Journal of biological engineering, 6:1–13, 2012.
  * [30] Joshua L Cherry and Frederick R Adler.  How to make a biological switch.  Journal of theoretical biology, 203(2):117–133, 2000.
  * [31] Jules Guilberteau, Camille Pouchol, and Nastassia Pouradier Duteil.  Monostability and bistability of biological switches.  Journal of Mathematical Biology, 83(6):65, 2021.
  * [32] Jie Li and Weinian Zhang.  Transition between monostability and bistability of a genetic toggle switch in escherichia coli.  Discrete & Continuous Dynamical Systems-Series B, 25(5), 2020.
  * [33] Simone Brabletz and Thomas Brabletz.  The zeb/mir-200 feedback loop—a motor of cellular plasticity in development and cancer?  EMBO reports, 11(9):670–677, 2010.
  * [34] Alexandra C Title, Sue-Jean Hong, Nuno D Pires, Lynn Hasenöhrl, Svenja Godbersen, Nadine Stokar-Regenscheit, David P Bartel, and Markus Stoffel.  Genetic dissection of the mir-200–zeb1 axis reveals its importance in tumor differentiation and invasion.  Nature communications, 9(1):4671, 2018.
  * [35] Menghan Chen and Ruiqi Wang.  Computational analysis of synergism in small networks with different logic.  Journal of Biological Physics, 49(1):1–27, 2023.
  * [36] Lawrence Perko.  Differential equations and dynamical systems.  2013\.
  * [37] Jinfang Zhu, Hidehiro Yamane, and William E Paul.  Differentiation of effector cd4 t cell populations.  Annual review of immunology, 28:445–489, 2009.
  * [38] Guiyuan Wang, Zhuoqin Yang, and Marc Turcotte.  Dynamic analysis of the time-delayed genetic regulatory network between two auto-regulated and mutually inhibitory genes.  Bulletin of Mathematical Biology, 82:1–30, 2020.
  * [39] Hongguang Xi and Marc Turcotte.  Parameter asymmetry and time-scale separation in core genetic commitment circuits.  Quantitative Biology, 3:19–45, 2015.
  * [40] Kiresh Parmar, Konstantin B Blyuss, Yuliya N Kyrychko, and Stephen John Hogan.  Time-delayed models of gene regulatory networks.  Computational and mathematical methods in medicine, 2015, 2015.
  * [41] Dongya Jia, Mohit Kumar Jolly, William Harrison, Marcelo Boareto, Eshel Ben-Jacob, and Herbert Levine.  Operating principles of tristable circuits regulating cellular differentiation.  Physical biology, 14(3):035007, 2017.
  * [42] Kishore Hari, Pradyumna Harlapur, Aditi Gopalan, Varun Ullanat, Atchuta Srinivas Duddu, and Mohit Kumar Jolly.  Emergent properties of coupled bistable switches.  Journal of Biosciences, 47(4):81, 2022.
  * [43] David S Glass, Xiaofan Jin, and Ingmar H Riedel-Kruse.  Nonlinear delay differential equations and their application to modeling biological network motifs.  Nature communications, 12(1):1788, 2021.
  * [44] Qing Hu, Min Luo, and Ruiqi Wang.  Identifying critical regulatory interactions in cell fate decision and transition by systematic perturbation analysis.  Journal of Theoretical Biology, 577:111673, 2024.
  * [45] Alexandra C Title, Sue-Jean Hong, Nuno D Pires, Lynn Hasenöhrl, Svenja Godbersen, Nadine Stokar-Regenscheit, David P Bartel, and Markus Stoffel.  Genetic dissection of the mir-200–zeb1 axis reveals its importance in tumor differentiation and invasion.  Nature communications, 9(1):4671, 2018.
  * [46] T Celià-Terrassa, C Bastian, DD Liu, B Ell, NM Aiello, Y Wei, J Zamalloa, AM Blanco, X Hang, D Kunisky, et al.  Hysteresis control of epithelial-mesenchymal transition dynamics conveys a distinct program with enhanced metastatic ability. nat commun. 2018; 9: 5005.