# Novel Generic Models for Differentiating Stem Cells
Reveal Oscillatory Mechanisms

Saeed Farjami  Department of Microbial Sciences, University of Surrey, GU2
7XH, Guildford, UK  Karen Camargo Sosa  Department of Biology and
Biochemistry, University of Bath, BA2 7AY, Bath, UK
Jonathan H.P. Dawes  Department of Mathematical Sciences, University of Bath,
BA2 7AY, Bath, UK  Robert N. Kelsh  Department of Biology and Biochemistry,
University of Bath, BA2 7AY, Bath, UK  Andrea Rocco 111Corresponding author:
<EMAIL_ADDRESS>Department of Microbial Sciences, University of Surrey,
GU2 7XH, Guildford, UK  Department of Physics, University of Surrey, GU2 7XH,
Guildford, UK

###### Abstract

Understanding cell fate selection remains a central challenge in developmental
biology. We present a class of simple yet biologically-motivated mathematical
models for cell differentiation that generically generate oscillations and
hence suggest alternatives to the standard framework based on Waddington’s
epigenetic landscape. The models allow us to suggest two generic dynamical
scenarios that describe the differentiation process. In the first scenario
gradual variation of a single control parameter is responsible for both
entering and exiting the oscillatory regime. In the second scenario two
control parameters vary: one responsible for entering, and the other for
exiting the oscillatory regime. We analyse the standard repressilator and four
variants of it and show the dynamical behaviours associated with each
scenario. We present a thorough analysis of the associated bifurcations and
argue that gene regulatory networks with these repressilator-like
characteristics are promising candidates to describe cell fate selection
through an oscillatory process.

Keywords: stem cell, multipotency, gene expression, gene regulatory network,
differentiation, oscillation

##  1 Introduction

In 1940, Waddington proposed representing the complex regulatory dynamics
driving the process of cellular differentiation as an ’epigenetic landscape’
[1] through which a single cell can be thought of as travelling. The
differentiating cell encounters successive ‘decision points’ in the landscape
morphology that correspond to differentiation events. These decision points
emerge from dynamical changes in the underlying cellular gene regulatory
network (GRN), and are mathematically described by variations in parameters
that influence the GRN dynamics. For example in zebrafish it is known that Wnt
signalling plays a fundamental role in specification and commitment of
melanocytes [2, 3, 4], and the level of Wnt expression may therefore be
considered as a parameter, driving the cellular GRN through one or more of the
decision points hypothesised by Waddington.

Despite the philosophical attractiveness of the epigenetic landscape metaphor,
the details have remained unclear and no completely self-consistent
mathematical description has emerged. For example there are debates regarding
the types of bifurcations characterising differentiation: saddle-node
bifurcations [5, 6] have been proposed as being more consistent with the
expected topologies of the core GRN responsible for differentiation than
pitchfork bifurcations. Experimental support for these mathematical
bifurcation phenomena, e.g. in embryonic stem cells (ESCs), is also lacking.
For example, when progesterone is washed away from mature Xenopus oocytes,
they do not dedifferentiate but remain mature [22], a finding that is more
consistent with the saddle-node bifurcation than with Waddington’s landscape.

Of course, the detailed biological understanding of cell differentiation
dynamics has grown significantly recently and it is appreciated that
qualitatively different mechanisms may be at work in different contexts. For
instance, detailed study of the differentiation of ESCs into two cell types
that form either embryonic or extraembryonic tissues indicates that this does
not occur in a simple deterministic manner [7]. Instead, cells appear to pass
transiently through intermediate states, each seemingly primed to
differentiate into one of the possible cell fates [8, 9]. Hence ESCs might
maintain a collection of meta-stable states, and a better framework to
understand and interpret experimental data is required [9].

Most fundamentally, the GRN topologies generating such dynamic stem cell
states are poorly understood. In the present work we focus on the question of
how broader multipotency (i.e. more than two fates) might generically be
generated, in a way that allows differentiation to be biased towards locally-
favourable outcomes. We note this is a topic of much current interest in
systems biology, particularly noting the exploration of connections between
GRN topology and oscillatory dynamics [12], generic links with differentiation
dynamics [13, 14, 15], and e.g. reshaping the epigenetic landscape by rewiring
the GRN using cross-repression among transcription factors (TFs) with self-
regulation [14].

In this paper, we analyse and compare five minimally-constructed GRNs that
admit cyclical dynamics and allow evolution towards fate commitment controlled
by an external stimulus. Our GRNs are variants of the standard repressilator,
a simple but fundamental circuit made up of three genes repressing each other.
Theoretical modelling of this circuit and its subsequent engineering in E.
coli have shown its capability to reproduce oscillatory behaviours under broad
parameter ranges [17].

We are here interested in adopting this model and the variants proposed to
mimick entry into the dynamical oscillatory state characterizing multipotency,
and the exit from it, towards a fully differentiated cell state.

Therefore, besides exhibiting oscillations, the system must be equipped with a
mechanism to exit the oscillatory regime and proceed towards differentiation.
Hence the system must exhibit a sequence of at least two bifurcations,
accounting for entering and exiting the oscillatory phase.

We identify two generic scenarios for entering and exiting the oscillatory
regime, and we test both in our proposed GRN models. In the first scenario
(S1), we hypothesise that a single parameter (i.e. one intervening signalling
pathway) drives the system both into and out of the oscillatory regime. An
example is Wnt signalling, required for the induction of neural crest, but
also for specification of both melanocyte and sensory neuron fates through
activation of key TFs (Mitf and Neurogenin, respectively) [3, 4, 18, 19, 20,
21, 22, 23], which we envisage as having this two-fold role. First, increasing
Wnt signalling promotes neural crest cells to enter an oscillatory multipotent
phase in which fate specification TFs are cyclically expressed; secondly, at
higher signalling levels, it promotes the cell’s exit from the oscillatory
phase allowing specification of single cell fates through increased Mitf or
Neurogenin expression which drive melanocyte or sensory neuron differentiation
respectively.

In the second scenario (S2) we hypothesise again that a first bifurcation to
oscillatory behaviour is driven by an external signalling pathway, but its
influence is then blocked above a critical value. Exit from the oscillatory
regime results from other signalling pathways, associated with other
parameters. This would arise for a multipotent progenitor that in response to
a specific signal starts oscillating, as in S1, driving expression of multiple
fate specification TFs, each characteristic of a different ’sub-state’. Close
to each transcriptional sub-state, the stem cell would express distinct cell
signalling receptors; then sufficient activation of these receptors by
environmental ligands could provide a mechanism that would move the cell out
of the oscillatory phase, thus driving differentiation.

The paper is structured as follows. In Section 2, we consider the simplest GRN
exhibiting oscillatory behaviour, the so-called repressilator [17]. We find
that the system exhibits a parameter regime where trajectories tend towards an
attracting limit cycle featuring a slowdown when passing near three sub-
states; this appears as a rather weak effect but is numerically detectable.
However, exit from the oscillatory regime leads to a single equilibrium,
making this GRN unable to select between alternative differentiated states.

This result motivates Sections 3 and 4 where we extend the standard
repressilator to include a second repressive circuit opposing the first; we
therefore term this GRN the ‘cross-repressilator’. Within this new GRN we
consider whether the two TFs act on each gene as either an ‘OR gate’ or an
‘AND gate’, leading to two variants of this circuit. We find that for the ‘OR
gate’, the circuit first transitions to oscillatory dynamics, as the standard
repressilator, but then followed by a further transition to a stably
differentiated cell state, in which only one of the genes is stably expressed.
The ‘AND gate’ is even more interesting, since it allows simultaneous
expression of two out of the three genes. In both cases our results appear to
hold for a large class of modelling choices for these gene interactions.

In Section 5 we extend the ‘AND gate’ cross-repressilator circuit to four and
five genes. We show that co-expression of up to three genes is possible
depending on the topology of the GRN. These dynamics are compatible with both
S1 and S2, as they provide an exit mechanism dependent on a single or multiple
bifurcation parameters into a number of alternative equilibria. Furthermore,
we note that in all the models, when a cell is in the oscillatory phase,
scoring for gene expression of fate-specific TFs would result (in a snapshot
view) in the cell being considered ‘fate specified’, despite the fact that it
actually retains full multipotency.

A discussion of the biological implications of our models is presented in
Section 6.

##  2 The standard repressilator

We first consider the repressilator circuit [17], whose GRN is shown in Fig.
1. Here each gene represents a master regulator transcription factor (TFs) of
the differentiation process for a specific fate.

Figure 1: A schematic representation of the repressilator circuit.

Assuming fast (un)binding dynamics of the TFs to DNA, and fast mRNA dynamics,
we can reduce the system by adiabatic elimination (see Supplementary
Materials) [24], and describe its dynamics with ordinary differential
equations (ODEs):

$\left\\{\begin{aligned}
\frac{dx_{1}}{dt}=b_{1}+\frac{g_{1}}{1+\alpha_{1}x_{2}^{h_{1}}}-d_{1}x_{1},\\\
\frac{dx_{2}}{dt}=b_{2}+\frac{g_{2}}{1+\alpha_{2}x_{3}^{h_{2}}}-d_{2}x_{2},\\\
\frac{dx_{3}}{dt}=b_{3}+\frac{g_{3}}{1+\alpha_{3}x_{1}^{h_{3}}}-d_{3}x_{3}.\end{aligned}\right.$
(1)

Here, $x_{i}$ ($i=1,\,2,\,3$) describe the concentrations of the proteins $\sf
X_{i}$, $b_{i}$ are background expression rates, $g_{i}$ are maximal gene
expression rates (which combine transcription, translation and mRNA
degradation rates), $\alpha_{i}$ are ratios of binding to unbinding rates
(association constants) of TFs to DNA, $h_{i}$ are Hill coefficients
describing possible cooperative effects (like multimerisation of the TFs
[24]), and $d_{i}$ are TF degradation rates.

Analysis (see e.g. [17]) shows that the repressilator supports oscillatory
behaviour over wide regions of parameter space. Fig. 2 shows (from left)
numerical solutions exhibiting the emergence of a family of stable periodic
orbits through a Hopf bifurcation as the parameter $g_{1}$ increases, time
series of sustained oscillations, and the limit cycle in phase space.
Oscillatory behaviour exists when $h>h^{*}\approx 2.20498$, for the parameter
values used in Fig. 2.

Figure 2: a) Bifurcation diagram for the repressilator as $g_{1}$ varies. A
family of stable periodic orbits (olive), shown with maxima and minima of
$x_{1}$, emanates from a supercritical Hopf bifurcation at $g_{a}$ and
terminates at another supercritical Hopf bifurcation at $g_{b}$. The stable
(solid) branch of equilibria (black) becomes unstable (dashed) between these
two bifurcations. The dotted vertical indicates the $g_{1}$-value used for
panels (b) and (c). b) Time courses of oscillatory concentrations $x_{1}$,
$x_{2}$, $x_{3}$ when $g_{1}=0.2$, $g_{2}=g_{3}=0.05$. c) The solution shown
in panel (b) in the phase space (grey) converging to the periodic orbit
(olive). The cross ($\times$) is the unique saddle equilibrium, at
$x_{eq}=(0.58469,\,0.87298,\,0.45578)$. Remaining parameter values are:
$b_{i}=10^{-5}$, $\alpha_{i}=50$, $d_{i}=0.01$, $h_{i}=3$, for $i=1,2,3$.

For $h=3$, the oscillations exist between the two Hopf bifurcations at
$g_{a}\approx 0.00764$ and $g_{b}\approx 15.7919$ suggesting that the system
is an example of S1 – variations in a single parameter are sufficient to enter
and also to exit the oscillatory regime. Our hypothesis for S1 is that
external signalling might increase $g_{1}=g_{1}(t)$ over time, letting the
system sequentially explore non-oscillatory, oscillatory, and a new non-
oscillatory regime. However, exit from the oscillatory regime at
$g_{1}(t)=g_{b}$ implies convergence onto the single available equilibrium,
making the GRN unable to capture the idea of multiple equilibria necessary to
describe stem cell differentiation.

Further, the period of the limit cycles remains bounded as $g_{1}$ varies. The
shape of the limit cycle (Fig. 2(c)) shows that the amplitudes of the
oscillations for different protein concentrations may differ significantly,
and generally appear to spike as the orbit approaches the coordinate axes; two
proteins are expressed at relatively low levels while the third is higher. Our
simulations show that even when the spiking of gene $\sf X_{1}$ is more
pronounced, with a huge increase in oscillation amplitude, only a logarithimic
or algebraic dependency of the oscillation period appears (see Supplementary
Materials), demonstrating the absence of any substantial slow-down around the
genes maximal expression substates, and failing thereby to effectively
describe any of the scenarios S1 or S2. Furthermore, no mechanism is provided
in this model for the stem cell to differentiate into more than one cell type,
as no alternative equilibria exist. Thus the standard repressilator does not
capture the qualitative changes in behaviour associated with differentiation.

##  3 The cross-repressilator with an ‘OR gate’

We now introduce a variant of the standard repressilator (‘cross-
repressilator’), where the inhibitory cycle is accompanied by a second
inhibitory cycle arranged in the opposite direction, as depicted in Fig. 3. We
assume that regulation of each gene follows an ‘OR gate’, i.e. each gene is
active when _both_ inhibitors are absent (see Supplementary Materials).

Figure 3: Illustration of cross-repression among three genes $\sf X_{1}$, $\sf
X_{2}$ and $\sf X_{3}$ formed with two copies of the circuit shown in Fig. 1,
oriented in opposing directions.

Making the same assumptions as the standard repressilator, the dynamics of
this GRN can be written as

$\left\\{\begin{aligned}
\frac{dx_{1}}{dt}=b_{1}+\frac{g_{1}}{(1+\alpha_{1}x_{2}^{h_{12}})(1+\beta_{1}x_{3}^{h_{13}})}-d_{1}x_{1},\\\
\frac{dx_{2}}{dt}=b_{2}+\frac{g_{2}}{(1+\alpha_{2}x_{3}^{h_{23}})(1+\beta_{2}x_{1}^{h_{21}})}-d_{2}x_{2},\\\
\frac{dx_{3}}{dt}=b_{3}+\frac{g_{3}}{(1+\alpha_{3}x_{1}^{h_{31}})(1+\beta_{3}x_{2}^{h_{32}})}-d_{3}x_{3},\end{aligned}\right.$
(2)

where the repression of each gene is described by the multiplication of two
decreasing Hill functions to capture the assumed ‘OR gate’, (see Supplementary
Materials). Like system Eq. 1, the last term in each equation represents
protein degradation. For simplicity, and since it does not cause qualitative
changes in the system dynamics, we remove basal expression here and from now
throughout this paper. Also, in the numerical simulations we set parameter
values equal for each gene, i.e. $g_{i}=g$, $\alpha_{i}=\alpha$,
$\beta_{i}=\beta$, $d_{i}=d$, and $h_{ij}=h$ for $i,j=1,2,3$; we return to
this in Section 3.2.

Numerical simulations of ODEs Eq. 2 are shown in Fig. 4. Figure 4(a) shows the
bifurcation diagram for system Eq. 2 with respect to the common parameter $g$,
all other parameters fixed. The (dashed) solid curves indicate (un)stable
equilibria (black) or (maxima and minima of) periodic orbits (olive green).
For $g\lessapprox 0.4$, system Eq. 2 has a single stable equilibrium which is
(due to the equalities in parameter values) symmetric: $x_{1}=x_{2}=x_{3}$.
The cyclic symmetry apparent in Eq. 2 implies that the Jacobian matrix
evaluated at this symmetric equilibrium has a pair of eigenvalues with equal
real parts (see Supplementary Materials). When this complex conjugate pair
crosses the imaginary axis, a supercritical Hopf bifurcation occurs generating
a family of stable periodic orbits. As $g$ increases further, the stable
periodic orbits disappear at a global bifurcation involving three new
equilibria created at saddle-node bifurcations.

Figure 4: a) Bifurcation diagram with respect to $g$ in log-scale for
$\alpha=9$, $\beta=0.1$, $h=3$, and $d=0.2$. b) The period of the periodic
orbits shown in panel (a) increases to infinity while $g$ tends to the SNIC
bifurcation. c) Time courses of the expression level of genes $\sf X_{1}$
(red), $\sf X_{2}$ (green) and $\sf X_{3}$ (blue) showing the sustained
oscillations, for $g=1.5$. d) For $g=1.6$, the system stops oscillating and
exits to a differentiated cell.

Moreover, there are two distinct global bifurcation mechanisms that remove the
periodic orbits. For $\alpha=9$ (shown in Fig. 4), we find that the periodic
orbits disappear at a SNIC (Saddle-Node on an Invariant Cycle) bifurcation
where three saddle-node bifurcations, related by the cyclic symmetry, occur on
the periodic orbit simultaneously. The SNIC bifurcation is well-known in
simple dynamical systems, such as the nonlinear pendulum with forcing and
damping [25] and in oscillators under external periodic perturbations [26,
27]. For $g>g_{\rm SNIC}\approx 1.51449$, system Eq. 2 has no periodic orbit
but three equilibria related by the cyclic symmetry. For each of the
equilibria, one of the $x_{i}$ is significantly larger than the other two.
Initial conditions (ICs) determine which of these equilibria attract the
trajectory.

Figure 4(b) shows that, in contrast with the standard repressilator, the
period of the periodic orbits diverges as $g$ approaches $g_{\rm SNIC}$.
Typical trajectories spend longer and longer in the vicinity of the three
‘slow regions’ of phase space where the saddle-node bifurcations are about to
occur, and exhibit rapid transitions between these ‘slow regions’. This
phenomenon has been referred to as the ‘ghost of SNIC’ [30], and indicates the
possibility of S2. We identify these slow regions with the sub-states
discussed earlier.

In Fig. 4 (right), we show time courses of system Eq. 2 before (c) and after
(d) the SNIC bifurcation with same ICs. Figure 4(c) shows the time evolution
of the three genes for g$g=1.5<g_{\rm SNIC}$. Whenever one of the genes is
highly expressed, the expression levels of the other two is low. In Fig. 4(d),
we observe the exit from the oscillatory behaviour associated with a further
increase of parameter $g$ ($g=1.6$). Here, gene $\sf X_{1}$ is stably
expressed at a high level, while genes $\sf X_{2}$ and $\sf X_{3}$ are stably
expressed at very low levels.

In summary, for $g$ increasing, the system shows three different behaviours:
(i) convergence to a (unique) stable equilibrium; (ii) oscillations among
three sub-states; and (iii) attraction to one of the stable equilibria. The
natural interpretation, compatible with S1, is to see the first two cases as
reflecting multipotency within the stem cell, while the third case represents
attainment of one of multiple differentiated cell types.

The phenotype attained in the differentiated state depends on the IC used for
the simulation and its intrinsic properties, here illustrated through
different $g$-values, as shown in Fig. 5. Figure 5(a) colour-codes the
$(x_{1},\,x_{2})$-plane according to the final equilibrium for three different
ICs for $x_{3}(t)$, and fixed $g$. Figure 5(b) shows the effect of varying
parameter $g$ on the same colourmap when $x_{3}(0)$ is fixed. The three
coloured regions in the $(x_{1},\,x_{2})$-plane twists around a central point.
Increasing $x_{3}(0)$ tightens the twists of these spiralling regions while
the opposite occurs when $g$ increases. Moreover, variations in $g$ keep the
central twisting point fixed at $(x_{1}(0),\,x_{2}(0))=(1,1)$ whereas
increasing $x_{3}(0)$ makes this point drift away.

Figure 5: a) Colourmaps indicating which gene is expressed starting from ICs
in the $(x_{1},\,x_{2})$-plane, for $g=2$ and three different ICs of $x_{3}$:
$x_{3}(0)=0.5,1$ and $1.5$. b) The same colourmaps as panel (a) for
$x_{3}(0)=1$ and three different values of $g$: $g=2,4$ and $6$.

To obtain insight into the feasibility of S2 for this system, we explore the
bifurcation behaviour in the $(g,\,\alpha)$-plane. Fig. 6(a) shows the
continuation of the bifurcations identified for $\alpha=9$ in the
$(g,\,\alpha)$-plane. When $\alpha$ decreases, the location of the Hopf
bifurcation curve ($\sf HB$) moves to larger $g$ whereas the curve of saddle-
node bifurcations ($\sf SN$) on which the SNIC bifurcation occurs moves to
lower values of g. At $\alpha=\alpha^{\ast}\approx 3$, the SNIC regime
terminates at a codimension-two bifurcation, indicated by a grey dot. Two new
bifurcation curves emerge from this point, indicated by the labels $\sf HC$
(homoclinic) and $\sf SN$ (saddle-node). For $\alpha<\alpha^{\ast}$ the three
saddle-node bifurcations (due to symmetry) occur away from the periodic orbit
in phase space, resulting in a region of the $(g,\,\alpha)$-plane in which
these new equilibria (three stable and three saddle equilibria) coexist with
the stable periodic orbits. The periodic orbits then collide with the three
saddle equilibria at a global bifurcation and disappear. As in the SNIC
bifurcation, the period of the periodic orbits goes to infinity as the orbits
approach the global bifurcation, but the mathematical details of the behaviour
differ from the SNIC case, as shown for example by the scaling of the period
with distance to the bifurcation point being different (see Supplementary
Materials).

Figure 6: a) Two-parameter continuation of saddle-node (black), Hopf (dark
green) and heteroclinic (red) bifurcations, shown in Fig. 4(a), in the
$(g,\alpha)$-plane; vertical axis ($\alpha$) is shown in log-scale. The dotted
and dashed lines indicate the corresponding $\alpha$ values used in Fig. 4(a)
($\alpha=9$) and in panel (b) ($\alpha=1$), respectively. b) Bifurcation
diagram of system Eq. 2 with respect to $g$ in log-scale when $\alpha=1$. The
diagram structure is mainly similar to the one in Fig. 4(a). However, here, a
small region of multi-stability is created between the oscillatory regime and
the stable equilibria such that the family of periodic orbits terminates at a
heteroclinic bifurcation (instead of a SNIC).

We observe also that as $\alpha$ decreases further the heteroclinic
bifurcation ($\sf HC$) and $\sf SN$ curves move further apart from each other
in the $(g,\alpha)$-plane, with maximum separation at $\alpha\approx 0.1$.

For $\alpha<\alpha^{\ast}$ the $\sf HB$ (green) and $\sf HC$ curves (red)
determine the left and right boundaries of the oscillatory region,
respectively. Also note that the $\sf HB$ and $\sf HC$ curves meet and are
tangent to each other at $\alpha=0.1$ which is not generic but arises here due
to the additional symmetry existing when $\alpha=0.1=\beta$. When
$\alpha=\beta$, ODEs Eq. 2 are symmetric under two independent symmetry
operations: the cyclic permutation symmetry
$(x_{1},x_{2},x_{3})\rightarrow(x_{2},x_{3},x_{1})$ and the interchange of any
two of the three variables, e.g.
$(x_{1},x_{2},x_{3})\rightarrow(x_{2},x_{1},x_{3})$. In group theory notation,
the symmetry group is now the symmetries of an equilateral triangle $D_{3}$
rather than just the cyclic group of order three, $\mathbb{Z}_{3}$. This
additional symmetry implies the tangency between the two curves. The analysis
around this point will be discussed in detail in a separate paper [32] and has
been observed in other contexts [28, 29].

Figure 6(b) shows the bifurcation diagram, varying $g$, for $\alpha=1$
corresponding to the lower horizontal dashed line indicated in Fig. 6(a). Here
we observe an interval of $g$-values in which the stable periodic orbit
coexists with the three stable equilibria. Consequently we would expect that
this transition to a stable equilibrium occurs rather smoothly since the sub-
states already exist prior to the termination of the oscillations.

These observations support S2. Once the system has entered the oscillatory
regime by tuning of $g$ through a first signalling pathway, a second
signalling pathway might intervene to lower $\alpha$, stop the oscillations,
and drive the system to differentiation. This scenario requires that the
change of $\alpha$ happens faster than the permanence of the system close to a
selected sub-state. This slowdown of the cycling dynamics is therefore
essential for this mechanism.

Thus, after establishment of the stem cell, increasing the first signal
initiates the cyclical expression of fate- specific TFs corresponding to the
different differentiated cell-types. Subsequently, in each primed sub-state,
as a result of TF specific activation of sensitivity to fate specification
signals, the cell is sensitised to local signals. When a cell receives such a
signal for sufficient time, this shifts $\alpha$ and removes the oscillations,
thereby initiating differentiation.

###  3.1 Illustrating the two scenarios with time-dependent parameters

So far we have assumed that all the parameters of system Eq. 2 are fixed and
do not vary in time. In this section, we explore more directly the two
scenarios S1 and S2 by assuming that relevant parameters are regulated in time
by an external signalling pathway. We will consider the case when either one
or both of $g$ and $\alpha$ are time dependent. In S1, parameter $g$ only is
responsible for shifting the cell from the multipotent regime into the
oscillatory regime and then into the differentiated cell. Accordingly, we
assume that $g(t)$ follows the form of an increasing Hill function, from 0 and
saturating at $\hat{g}>g_{\rm SNIC}$ as

$g(t)=\frac{\hat{g}t^{p}}{\tau^{p}+t^{p}},$

where we set $\hat{g}=1.7$, $p=1.25$ and $\tau=1000$. The bottom panel of Fig.
7(a) illustrates the time-variation of $g(t)$. The top panel in Fig. 7(a)
shows how the system responds to $g(t)$, moving through the three different
regimes as $g(t)$ increases and tends to the value $\hat{g}$ (compare with
Fig. 4(a)). In detail, the system first settles quickly from the IC to the
single stable equilibrium and follows it as $g(t)$ slowly approaches the Hopf
bifurcation. After Hopf we observe the onset of rapid but finite-frequency
oscillations. As $g\longrightarrow g_{\rm SNIC}$, the period of oscillations
increases significantly, giving rise to dynamical phases close to each of the
three sub-states. Finally, when $g(t)$ crosses $g_{\rm SNIC}$, the
oscillations disappear and the system settles at one of the three equilibria,
each corresponding to a different expressed gene. For example, in Fig. 7, gene
$\sf X_{3}$ (blue) is expressed and the other two genes remain at a low level.

Figure 7: a) The development over time from a multipotent cell to a
differentiated cell controlled by a time-dependent profile $g(t)$ (shown in
the lower panel) according to S1. Parameters are as Fig. 4(c). b) Development
from a stable oscillatory regime into a differentiated cell according to S2,
when the secondary parameter $\alpha$ starts decreasing at $\tilde{t}=0$
(shown in the lower panel). Here $g=1.5$ is held constant, and all other
parameters take the same values as in Fig. 4(c).

As discussed for S1, the precise state of the system when $g$ goes through
$g_{\rm SNIC}$, as well as the value of the saturation level $\hat{g}$ play
crucial roles in selecting cell fate. The shape of the profile of $g(t)$ is
also relevant for this selection.

For S2, we hypothesise that two parameters, $g$ and $\alpha$, are involved in
cell differentiation. As shown in Fig. 7(b), we first assume that parameter
$g$ shifts the system into the oscillatory regime while $\alpha$ remains
fixed. At a certain time, we suppose a second pathway causes a decrease in
$\alpha$, modelled as

$\alpha(t)=9-\mathcal{H}(t-\tilde{t})\times\frac{(t-\tilde{t})^{p}}{\tau^{p}+(t-\tilde{t})^{p}},$

where $\mathcal{H}$ is a Heaviside step function (i.e. $\mathcal{H}(s)=0$ when
$s<0$ and $\mathcal{H}(s)=1$ when $s>0$), $\tilde{t}$ is the time when the
time dependency of $\alpha(t)$ starts, and the constants $\tau=1000$ and $p=1$
control the profile shape.

S2 clearly allows more control over the final fate of the cell through the
increased number of coefficients involved in specifying the time-dependency of
the second bifurcation parameter. To understand these further dependencies, we
analyse the effect of variations in the switch-on time $\tilde{t}$ and the
time scale parameter $\tau$ for the saturation of $\alpha(t)$. Shown in Fig.
8, we colour-code the parameter plane $(\tau,\tilde{t})$ based on the
expressed gene for $g=1.5$ (a) and $g=1.51$ (b). The $(\tau,\tilde{t})$-plane
is divided into stripes which periodically alternate between red, green and
blue. The lower value of $g$ in Fig. 8(a) corresponds to narrower stripes as a
result of the smaller oscillation period, with the system being more sensitive
to changes in switch-on time and saturation time scale. Further, more
horizontal stripes in Fig. 8(b) indicate that cell fates become less sensitive
to $\tau$ for this larger value of $g$, at least in the regime close to the
global bifurcation.

Figure 8: Colourmaps show the expressed gene in $(\tau,\tilde{t})$-plane
starting from the same IC for $g=1.5$ (a) and $g=1.51$ (b).

###  3.2 Breaking the symmetry

Until now, we assumed that system Eq. 2 is symmetric and all the association
rates $\alpha_{i}$ and $\beta_{j}$, $i,j=1,2,3$ for the TFs $\sf X_{1}$, $\sf
X_{2}$ and $\sf X_{3}$ are equal. Here, we break the symmetry in the
reciprocal repressing loops in Fig. 3 and consider the case that $\alpha_{i}$
and $\beta_{j}$ take different values, as specified in Fig. 9. Note that we
keep the remaining parameters $g$ and $d$, unchanged.

Figure 9(a) shows the bifurcation diagram with respect to $g$ in log-scale.
The structure remains qualitatively unchanged compared to the symmetric case.
For $g$ small, a unique stable equilibrium exists, which undergoes a Hopf
bifurcation and becomes unstable. For $g$ large, there are six branches of
equilibria; the stable branches merge with the saddle-type ones and disappear
as $g$ decreases. The family of stable periodic orbits emerging from the Hopf
bifurcation terminates at a SNIC bifurcation. However, here, only one
equilibrium lies on the invariant cycle contrary to the symmetric case, with
three equilibria on the invariant cycle. This happens because the saddle-node
bifurcations, at which the six branches of equilibria appear, occur at
different $g$-values (they are no longer symmetry-related), indicated by the
vertical dotted lines.

Figure 9: a) Bifurcation diagram of the system for $\alpha_{1}=7$,
$\alpha_{2}=9$, $\alpha_{3}=8$, $\beta_{1}=0.115$, $\beta_{2}=0.07$,
$\beta_{3}=0.1$ and $d=0.2$ in log-scale. b) Time courses of the oscillations
in the expression level of genes $\sf X_{1}$ (red), $\sf X_{2}$ (green) and
$\sf X_{3}$ (blue) for the same parameters in panel (a) and $g=1.4$.

Figure 9(b) shows time courses of the genes concentrations just before the
SNIC bifurcation at $g_{\rm SNIC}\approx 1.41054$; the other parameter values
are as in Fig. 9(a). As before, each TF is expressed cyclically when the other
two are at low levels. However, gene $\sf X_{2}$ remains expressed for a
significantly longer time than the other two, particularly $\sf X_{3}$. Due to
the occurrence of the saddle-node bifurcation with high level of $\sf X_{2}$
at smaller values of $g$, the slow region created around the bifurcation has a
stronger effect which leads to a longer ‘dwelling time’ at expressed $\sf
X_{2}$. In fact, the dwelling time for each expressed gene depends on the
distance of g-value at which the corresponding saddle-node bifurcation
(indicated with a dotted vertical line of the same colour) occurs.

Analysis of Fig. 9 determines that with the current parameters, the conditions
favour gene $\sf X_{2}$; therefore, when $g$ crosses to the right of the SNIC
bifurcation, the system settles at a stable equilibrium where $\sf X_{2}$ is
expressed. Changing this condition results in a shift in the location of the
saddle-node bifurcations along the horizontal axis and might alter the order
of saddle-node bifurcations; consequently, lengthening the dwelling time for
another gene.

##  4 The cross-repressilator with an ‘AND gate’

The third model we consider here is similar to system Eq. 2 with the
difference that each gene is repressed only if all other genes are highly
expressed (a so-called ‘AND gate’).

This model is described by

$\left\\{\begin{aligned}
\frac{dx_{1}}{dt}=b_{1}+\frac{g_{11}+g_{12}x_{2}^{h_{12}}+g_{13}x_{3}^{h_{13}}}{(1+\alpha_{1}x_{2}^{h_{12}})(1+\beta_{1}x_{3}^{h_{13}})}-d_{1}x_{1},\\\
\frac{dx_{2}}{dt}=b_{2}+\frac{g_{21}+g_{22}x_{3}^{h_{23}}+g_{23}x_{1}^{h_{21}}}{(1+\alpha_{2}x_{3}^{h_{23}})(1+\beta_{2}x_{1}^{h_{21}})}-d_{2}x_{2},\\\
\frac{dx_{3}}{dt}=b_{3}+\frac{g_{31}+g_{32}x_{1}^{h_{31}}+g_{33}x_{2}^{h_{32}}}{(1+\alpha_{3}x_{1}^{h_{31}})(1+\beta_{3}x_{2}^{h_{32}})}-d_{3}x_{3},\end{aligned}\right.$
(3)

where again the decreasing Hill functions describe the effect of repressing
genes (see Supplementary Materials). We also set $g_{ij}=g_{j}$ for
$i,j=1,\ldots,3$ and $d_{1}=d_{2}=d_{3}=d$. We also fix
$\alpha_{i}=\beta_{i}=1$ and the Hill function exponents $h_{ij}=3$ for
$i,j=1,2,3$.

As previously, we explore system Eq. 3 numerically and compute the bifurcation
diagram, varying $g_{1}$ while keeping $g_{2}$ and $g_{3}$ fixed (and
unequal). Figure 10 (a) shows a bifurcation diagram for system Eq. 3 with
respect to $g_{1}$. For small values of $g_{1}$, the system has a single
stable equilibrium (black solid line) which undergoes a supercritical Hopf
bifurcation and becomes unstable (dashed line). The stable periodic orbits
(olive) are indicated by the maximum and minimum values of $x_{1}$. The
periodic orbit, similar to the ‘OR gate’ case, terminates at a SNIC
bifurcation where it collides with a saddle-node bifurcation. Note that these
branches are part of a single branch of equilibria terminating at another
oscillatory regime at very large parameter values $g_{1}$ (not shown). Panels
(b) and (c) in Fig. 10 show time courses of the expression level of each gene
in system Eq. 3 for values of $g_{1}$ mentioned in the caption.

Figure 10: : a) Bifurcation diagram of system Eq. 4 with respect to $g_{1}$
for $g_{2}=0.76$, $g_{3}=1.9$ and $d=0.3$. b) Time courses of the expression
level of genes $\sf X_{1}$ (red), $\sf X_{2}$ (green) and $\sf X_{3}$ (blue)
showing the emergence of an intermediate expression level in the oscillatory
regime of system Eq. 3 when $g_{1}=2.8$. c) Gene expression levels converging
to equilibrium in the tristable regime for $g_{1}=3$.

Clearly, the behaviour of the ‘AND gate’ model follows closely the analysis
presented for system Eq. 2, but with one important and biologically relevant
difference. Compared to system Eq. 2, in which only one gene could be highly
expressed and the other two genes were necessarily at very low levels, in
system (3), two genes can be expressed at high or moderate levels
simultaneously. Of course, the precise levels depend on parameter values, but
the ‘AND gate’ structure is crucial in admitting states with more than one
high expression level. The form of the periodic orbit also changes and shows
that it lingers in sub-states where two gene expression levels are non-zero.
This notable feature is potentially useful in matching models with
quantitative experimental data where it allows one to select a model that
allows or prevents the existence of phenotypes dependent upon the expression
of a single or a combination of master regulator TFs.

##  5 Extension to four and five genes

So far, we considered GRNs with only three genes. We now investigate the
dynamics with more than three genes. For the same cross-repressilator in ‘AND
gate’ configuration discussed in the previous section, an odd number of genes
is necessary to create a negative feedback loop that supports sustained
oscillations [31]. Hence a trivial extension of the cross-repressilator to
four genes, preserving the cyclical topology, fails to reproduce oscillatory
dynamics, and therefore we do not consider it here. Even numbers of genes in
the cycle result in states where the gene expression levels alternate around
the loop. However, adding cross connections to this four-gene network so as to
build three-gene negative loops, embedded within the four-gene GRN as shown in
Fig. 11(a), recovers the oscillatory dynamics. For five genes, the simple
cyclical topology of the cross-repressilator, with only ‘nearest neighbour’
interactions and without additional cross-regulation, shown in Fig. 11(b),
maintains oscillatory behaviour.

Figure 11: a) Extension of the cross-repressilator GRN to four genes with
cross connections that support oscillations. b) A cyclical arrangement of five
genes, sufficient to drive oscillations without additional cross connections.

###  5.1 The four gene network

Here, we describe the dynamics of a GRN with four genes, as illustrated in
Fig. 11(a). The network is fully connected such that each gene is repressed by
the other three genes, operating in ‘AND gate’ configuration. The results in
(see Supplementary Materials for details):

$\left\\{\begin{aligned}
\frac{dx_{1}}{dt}=b_{1}+\frac{g_{11}+g_{12}x_{2}^{h_{12}}+g_{13}x_{3}^{h_{13}}+g_{14}x_{4}^{h_{14}}+g_{15}x_{2}^{h_{12}}x_{3}^{h_{13}}+g_{16}x_{2}^{h_{12}}x_{4}^{h_{14}}+g_{17}x_{3}^{h_{13}}x_{4}^{h_{14}}}{(1+\alpha_{1}x_{2}^{h_{12}})(1+\beta_{1}x_{3}^{h_{13}})(1+\gamma_{1}x_{4}^{h_{14}})}-d_{1}x_{2},\\\
\frac{dx_{2}}{dt}=b_{2}+\frac{g_{21}+g_{22}x_{3}^{h_{23}}+g_{23}x_{4}^{h_{24}}+g_{24}x_{1}^{h_{21}}+g_{25}x_{3}^{h_{23}}x_{4}^{h_{24}}+g_{26}x_{3}^{h_{23}}x_{1}^{h_{21}}+g_{27}x_{4}^{h_{24}}x_{1}^{h_{21}}}{(1+\alpha_{2}x_{3}^{h_{23}})(1+\beta_{2}x_{4}^{h_{24}})(1+\gamma_{2}x_{1}^{h_{21}})}-d_{2}x_{2},\\\
\frac{dx_{3}}{dt}=b_{3}+\frac{g_{31}+g_{32}x_{4}^{h_{34}}+g_{33}x_{1}^{h_{31}}+g_{34}x_{2}^{h_{32}}+g_{35}x_{4}^{h_{34}}x_{1}^{h_{31}}+g_{36}x_{4}^{h_{34}}x_{2}^{h_{32}}+g_{37}x_{1}^{h_{31}}x_{2}^{h_{32}}}{(1+\alpha_{3}x_{4}^{h_{34}})(1+\beta_{3}x_{1}^{h_{31}})(1+\gamma_{3}x_{2}^{h_{32}})}-d_{3}x_{3},\\\
\frac{dx_{4}}{dt}=b_{4}+\frac{g_{41}+g_{42}x_{1}^{h_{41}}+g_{43}x_{2}^{h_{42}}+g_{44}x_{3}^{h_{43}}+g_{45}x_{1}^{h_{41}}x_{2}^{h_{42}}+g_{46}x_{1}^{h_{41}}x_{3}^{h_{43}}+g_{47}x_{2}^{h_{42}}x_{3}^{h_{43}}}{(1+\alpha_{4}x_{1}^{h_{41}})(1+\beta_{4}x_{2}^{h_{42}})(1+\gamma_{4}x_{3}^{h_{43}})}-d_{4}x_{4}.\end{aligned}\right.$
(4)

For simplicity, we again set $g_{ij}=g_{j}$ and $d_{i}=d$ for
$i,j=1,\ldots,4$. We also keep fixed the values $g_{j}=0.5$ for
$j=5,\ldots,7$, $\alpha_{i}=\beta_{i}=\gamma_{i}=1$ and $h_{ij}=3$ for
$i,j=1,\ldots,4$ ($i\neq j$).

Figure 12(a) shows the bifurcation diagram of system Eq. 4 with respect to
$g_{1}$, together with time courses of the expression levels of all genes in
the system for two specific values of parameter $g_{1}$ in the multipotent (b)
and differentiated (c) regimes. In Fig. 12(a), a single branch of stable
equilibria (solid; black) undergoes a Hopf bifurcation and becomes
unstable/saddle (dashed). The envelope of stable periodic orbits (olive)
emanating from the Hopf bifurcation terminates at a SNIC bifurcation where new
branches of equilibria appear. Four of these new branches are stable, and
collide with the other four unstable branches, disappearing at saddle-node
bifurcations.

Figure 12: a) Bifurcation diagram of system Eq. 4 with respect to $g_{1}$ for
$g_{2}=1.6$, $g_{3}=0.5$, $g_{4}=1.5$ and $d=0.4$. (b) and (c) Time courses of
the expression level of genes $\sf X_{1}$ (red), $\sf X_{2}$ (green), $\sf
X_{3}$ (blue) and $\sf X_{4}$ (magenta) with two different behaviours: in the
multipotent oscillatory regime of system Eq. 4 when $g_{1}=3$ (b) and in the
multistable regime for $g_{1}=3.6$ (c).

Panels (b) and (c) in Fig. 12 show time series of the expression levels of the
genes in the oscillatory and differentiated regimes for $g_{1}=3$ and
$g_{1}=3.6$, respectively. One notable, and biologically relevant, feature of
this model is that all four genes can be expressed, with one gene highly
expressed, two genes expressed at intermediate levels, and one gene only
slightly expressed. This scenario is consistent with differentiation
reflecting a combinatorial role for multiple TFs. The relative concentrations
of the expressed genes would open up possibilities for differentiation to
multiple cell types.

###  5.2 The five gene network

We now consider the five-gene network depicted in Fig. 11(b). The network
comprises two inhibitory loops arranged in opposing directions so that gene
$X_{i}$ is repressed by genes $X_{i-1}$ and $X_{i+1}$, again in ‘AND gate’
configuration. Therefore, similar to the three-gene case, the dynamical
equations are

$\left\\{\begin{aligned}
\frac{dx_{1}}{dt}=b_{1}+\frac{g_{11}+g_{12}x_{2}^{h_{12}}+g_{13}x_{5}^{h_{15}}}{(1+\alpha_{1}x_{2}^{h_{12}})(1+\beta_{1}x_{5}^{h_{15}})}-d_{1}x_{1},\\\
\frac{dx_{2}}{dt}=b_{2}+\frac{g_{21}+g_{22}x_{3}^{h_{23}}+g_{23}x_{1}^{h_{21}}}{(1+\alpha_{2}x_{3}^{h_{23}})(1+\beta_{2}x_{1}^{h_{21}})}-d_{2}x_{2},\\\
\frac{dx_{3}}{dt}=b_{3}+\frac{g_{31}+g_{32}x_{4}^{h_{31}}+g_{33}x_{2}^{h_{32}}}{(1+\alpha_{3}x_{4}^{h_{34}})(1+\beta_{3}x_{2}^{h_{32}})}-d_{3}x_{3},\\\
\frac{dx_{4}}{dt}=b_{4}+\frac{g_{41}+g_{42}x_{5}^{h_{43}}+g_{43}x_{3}^{h_{43}}}{(1+\alpha_{4}x_{5}^{h_{45}})(1+\beta_{4}x_{3}^{h_{43}})}-d_{4}x_{4},\\\
\frac{dx_{5}}{dt}=b_{5}+\frac{g_{51}+g_{52}x_{1}^{h_{51}}+g_{53}x_{4}^{h_{54}}}{(1+\alpha_{5}x_{1}^{h_{51}})(1+\beta_{5}x_{4}^{h_{54}})}-d_{5}x_{5}.\end{aligned}\right.$
(5)

where the maximal expression rates imposed by the anticlockwise and the
clockwise loops respectively equal to $g_{1}$, $g_{2}$ and $d_{i}$ to $d$. We
also fix the values $\alpha_{i}=\beta_{i}=\gamma_{i}=1$ and $h_{ij}=3$ for
$i,j=1,\ldots,5$ ($i\neq j$). Figure 13(a) shows the bifurcation diagram
remains unchanged with respect to $g_{1}$ except for the number of stable
(solid; black) and saddle (dashed; black) branches in the multistable regime.
Figure 13(b) exhibits time courses of gene expression levels when $g_{1}=2$.
Figure 13(c) shows gene concentrations for $g_{1}=2.5$ when the cell
differentiates. In the differentiated regime, the number of expressed genes
increases to three compared with the three-gene case. We conjecture that more
generally for a cyclic network of $n$ genes, $\frac{n+1}{2}$ genes are
expressed at high or moderate (but unequal) levels and $\frac{n-1}{2}$ are
either not expressed at all or are only weakly expressed.

Figure 13: (a) Bifurcation diagram of system (5) with respect to $g_{1}$ when
$g_{2}=0.35$, $g_{3}=0.2$ and $d=0.4$. (b) Time series of the solutions in the
oscillatory/multipotent regime when $g_{1}=2$. (c), and multistable/
differentiated regime for $g_{1}=2.5$.

Biologically, the stable generation of cells expressing mixtures of key fate
specification TFs shown by these four- and five-gene models matches closely
the concept of combinatorial TF expression defining similar cell-types. This
is most keenly demonstrated in the peripheral and central nervous systems
where distinct neuron types express mixtures of TFs that together specify
their individual fate [33]. It may be that the fully differentiated phenotype
represents the impact of a mix of TFs expressed at different levels, even
where we tend to think of a single dominant master regulator, such as Mitf in
melanocyte development.

##  6 Conclusions

We have presented a detailed bifurcation analysis of four variants of the
standard repressilator model, and discussed their suitability as novel generic
models for the process of cellular differentiation. While the original
repressilator fails to provide convincing dynamical behaviour, we find that
the variants proposed here exhibit biologically relevant features. Our
modelling moves beyond Waddington’s vision of the epigenetic landscape since
such a picture cannot capture highly dynamical internal states of cells which
our model places centre stage.

Bifurcation analysis of the systems using [16] reveals that the oscillatory
behaviour emerges from a Hopf bifurcation. At the late stages of multipotency,
the oscillating cell visits a set of sub-states where the relevant fate-
specific transcription factor(s) are highly expressed relative to the other
gene(s). The oscillations terminate at a global bifurcation and the cell
settles at a state resembling that of a committed cell-type. From a biological
perspective, we can think of the individual fate-specific transcription
factors (TFs) as defining individual cell-fates. Where they oscillate, the
stem cell retains multipotency (it sequentially transits nearby all states),
but at the same time it becomes periodically biased to adopt individual fates
due, we propose, to changing exposure to fate-specification signals. When such
a signal is received at a sufficiently high level, it may force oscillations
to cease and a single TF becomes stably expressed at the expense of the others
–– cell differentiation has begun.

We have focused on two possible scenarios for driving the oscillatory
behaviour. In the first scenario, the same pathway is responsible for
initiating and terminating the oscillations, by the increase of the chosen
parameter. In the second scenario, a first pathway brings the system in the
oscillatory region, and a second one intervenes to drive the exit from it. Our
results show that the second scenario, however, gives more control over cell
fate. Moreover, the timing of controlling signals is essential in our
framework, since this combines with ICs to determine the specific TF that
remains stably expressed, and thus the differentiated cell-type.

Depending on the type of gates assumed in the network, one or two genes (or
more in networks of higher dimensionality) can be co-expressed during both the
oscillatory phase and in the selected differentiated cell type. During the
oscillations, combinations of genes may be expressed transiently –– these
might correspond to partially-restricted cell intermediates (expressing a
subset of fate-specific TFs); however, we emphasise that in our scenarios
these retain oscillations and hence must be thought of as fully multipotent,
despite their ‘snapshot appearance’. In the selected differentiated cell-type,
the combination and the levels of expressed TFs vary with ICs and timing of
the fate-specification factor; thus different cell-types, with specific
quantitatively distinguishable combinations of TFs, are generated. Such a
model now provides a view of how multipotency might be exhibited in stem
cells, in contrast to the series of bipotential intermediates that underpin
standard thinking in development. Note that the external signals drive a fully
multipotent cell to a specific committed state directly, without further
intermediates with restricted potency. Finally, the case with which stable
(committed) states consisting of mixed expression of two or more TFs can be
generated may well relate to the mechanism producing related cell-types, e.g.
neuron types in both the central and peripheral nervous system.

There is a wider context for the results presented here; the identification of
dynamical scenarios typical for repressilator-type GRNs, particularly in view
of their proposed role in clock circuits underlying global patterning
processes in development [34]. Although the developmental context considered
there is different, the formulation and behaviour of the mathematical models
has some points of similarity with our analysis, such as the occurrence of
Hopf and SNIC bifurcations, which are generic bifurcations involving time-
periodic oscillations.

A key difference with our work is ordering of the different regimes. Model 1
of [34] shows a progression from oscillations to a single stable equilibrium,
and then a saddle-node bifurcation that creates multiple equilibria. Model 2
of [34] has no Hopf bifurcation but instead a SNIC separating oscillatory from
multistable regime; there is no single equilibrium. In contrast, our model
robustly indicates a progression from a single equilibrium to oscillations via
a Hopf bifurcation, with the oscillations then terminating at a global
bifurcation where multiple equilibria emerge. We note that the distinction
between the SNIC mechanism and the combination of the global and saddle-node
bifurcations is likely to be extremely difficult to detect experimentally due
to the high level of resolution required in time, and the effects of noise. It
is however mathematically important in terms of comparisons between models.

In terms of model construction a key point of difference is that both Model 1
and Model 2 in [34] are constructed by interpolating, rather artificially,
between two mechanisms potentially driven by distinct and unrelated
morphogens. Our dynamical model is structurally more robust and parsimonious
as our work has shown that the transition from equilibrium to oscillations
whose frequency drops to zero, to ‘differentiated’ equilibria, is a natural
result of a single GRN mechanism. Naturally, spatial patterning is not part of
our model since we focus on the dynamics within a single cell, but
establishing connections between temporal dynamics and spatial patterning is
an intriguing research direction that we will explore in a future paper.

## Acknowledgments

We gratefully acknowledge Laura Hattam and Michael Thomas for helpful
preliminary modelling work. RNK gratefully acknowledges support from the
Institute for Mathematical Innovation at the University of Bath that helped to
advance this project in its early stages. Funding was provided by BBSRC grants
BB/S01604X/1 (SF and AR) and BB/S015906/1 (KCS, JHPD, RNK).

## References

  * [1] Waddington CH. 1940 _Organisers and Genes_. Cambridge, UK: Cambridge University Press.
  * [2] Vibert L, Aquino G, Gehring I, Subkhankulova T, Schilling TF, Rocco A, Kelsh RN. 2016 An ongoing role for Wnt signaling in differentiating melanocytes in vivo. Pigment Cell Melanoma Res. 30, 219–232. (doi:10.1111/pcmr.12568).
  * [3] Dorsky RI, Raible DW, Moon RT. 2000 Direct regulation of nacre, a zebrafish MITF homolog required for pigment cell formation, by the Wnt pathway. Genes & Dev. 14, 158–162. (doi:10.1101/gad.14.2.158).
  * [4] Lewis JL, Bonner J, Modrell M, Ragland JW, Moon RT, Dorsky RI, Raible DW. 2004 Reiterated Wnt signaling during zebrafish neural crest development. Development 131, 1299–1308. (doi:10.1242/dev.01007).
  * [5] Ferrell Jr JE. 2012 Bistability, bifurcations, and Waddington’s epigenetic landscape. Curr. Biol. 22, R458. (doi:10.1016/j.cub.2012.03.045).
  * [6] Huang S, Li F, Zhou JX, Qian H. 2017 Processes on the emergent landscapes of biochemical reaction networks and heterogeneous cell population dynamics: differentiation in living matters. J. R. Soc. Interface 14, 20170097. (doi:10.1098/rsif.2017.0097).
  * [7] Huang S. 2010 Cell Lineage Determination in State Space: A Systems View Brings Flexibility to Dogmatic Canonical Rules. PloS Biol. 8, e1000380. (doi:10.1371/journal.pbio.1000380).
  * [8] Wigger M, Kisielewska K, Filimonow K, Plusa B, Maleszewski M, Suwińska A. 2017 Plasticity of the inner cell mass in mouse blastocyst is restricted by the activity of FGF/MAPK pathway. Sci. Rep. 7, 15136 (doi:10.1038/s41598-017-15427-0).
  * [9] Graf T, Stadtfeld M. 2008 Heterogeneity of embryonic and adult stem cells. Cell Stem Cell 3, 480–483 (doi:10.1016/j.stem.2008.10.007).
  * [10] Nikaido M. et al., 2021 Zebrafish pigment cells develop directly from highly multipotent progenitors. Submitted.
  * [11] Kelsh RN, Camargo Sosa K, Farjami S, Makeev V, Dawes JHP, Rocco A. 2021 Cyclical Fate Restriction, a new view of neural crest cell fate specification. Submitted.
  * [12] Perez-Carrasco R, Barnes CP, Schaerli Y, Isalan M, Briscoe J, Page KM. 2018 Combining a Toggle Switch and a Repressilator within the AC-DC Circuit Generates Distinct Dynamical Behaviors. Cell Syst. 6, 521-530.e3. (doi:10.1016/j.cels.2018.02.008).
  * [13] Suzuki N, Furusawa C, Kaneko K, 2011 Oscillatory Protein Expression Dynamics Endows Stem Cells with Robust Differentiation Potential. Plos One 6, e27232. (doi:10.1371/journal.pone.0027232).
  * [14] Rabajante JF, Babierra AL. 2015 Branching and oscillations in the epigenetic landscape of cell-fate determination. Prog. Biophys. Mol. Biol. 117, 240–249. (doi:10.1016/j.pbiomolbio.2015.01.006).
  * [15] Kaity B, Sarkar R, Chakrabarti B, Mitra MK. 2018 Reprogramming, oscillations and transdifferentiation in epigenetic landscapes. Sci. Rep. 8, 7358 (doi:10.1038/s41598-018-25556-9).
  * [16] Doedel EJ, Champneys AR, Fairgrieve T, Kuznetsov YA, Oldeman BE, Paffenroth R, Sandstede B, Wang X, Zhang C. 2010 Auto-07p: Continuation and Bifurcation Software for Ordinary Differential Equations. Concordia University, Canada. Available from: http://cmvl.cs.concordia.ca/auto.
  * [17] Elowitz MB, Leibler S. 2000 A synthetic oscillatory network of transcriptional regulators. Nature 403, 335–338. (doi:10.1038/35002125).
  * [18] Saint-Jeannet JP, He X, Varmus HE, Dawid IB. 1997 Regulation of dorsal fate in the neuraxis by Wnt-1 and Wnt-3a. PNAS 94, 13713–13718. (doi:10.1073/pnas.94.25.13713).
  * [19] Dorsky RI, Moon RT, Raible DW. 1998 Control of neural crest cell fate by the Wnt signalling pathway. Nature 396, 370–373. (doi:10.1038/24620).
  * [20] Dunn KJ, Williams BO, Li Y, Pavan WJ. 2000 Neural crest-directed gene transfer demonstrates Wnt1 role in melanocyte expansion and differentiation during mouse development. PNAS 97, 10050–10055. (doi.org/10.1073/pnas.97.18.10050).
  * [21] Takada R, Takada S, Takeda K, Yasumoto K, Watanabe K, Udono T, Saito H, Takahashi K, Shibahara S. 2000 Induction of Melanocyte-specific Microphthalmia-associated Transcription Factor by Wnt-3a. The Journal of Biological Chemistry 275, 14013–14016. (doi:10.1074/jbc.C000113200).
  * [22] Jin EJ, Erickson CA, Takada S, Burrus LW. 2001 Wnt and BMP signaling govern lineage segregation of melanocytes in the avian embryo. Dev. Biol. 233, 22–37. (doi:10.1006/dbio.2001.0222).
  * [23] Lee HY, Kléber M, Hari L, Brault V, Suter U, Taketo MM, Kemler R, Sommer L. 2004 Instructive Role of Wnt/$\beta$-Catenin in Sensory Fate Specification in Neural Crest Stem Cells. Science 303, 1020–1023. (doi:10.1126/science.1091611).
  * [24] Greenhill ER, Rocco A, Vibert L, Nikaido M, Kelsh RN. 2011 An Iterative Genetic and Dynamical Modelling Approach Identifies Novel Features of the Gene Regulatory Network Underlying Melanocyte Development. PLoS Genet. 7, e1002265. (doi:10.1371/journal.pgen.1002265).
  * [25] Coullet P, Gilli JM, Monticelli N, Vandenberghe N. 2005 A damped pendulum forced with a constant torque. Am. J. Phys., 73, 1122–1128. (doi:10.1119/1.2074027).
  * [26] Andronov AA, Vitt AA, Khaikin SE. 1996 _Theory of Oscillators_. New York, US: Pergamon.
  * [27] Tsai TL, Dawes JHP. 2013 Dynamics near a periodically-perturbed robust heteroclinic cycle. Physica D 262, 14–34. (doi:10.1016/j.physd.2013.07.009).
  * [28] Gambaudo JM. 1985 Perturbation of a Hopf bifurcation by an external time-periodic forcing. J. Diff. Eq. 57, 172–199. (doi:10.1016/0022-0396(85)90076-2).
  * [29] Marques F, Meseguer A, Lopez JM, Pacheco JR, Lopez JM. 2013 Bifurcations with imperfect SO(2) symmetry and pinning of rotating waves. Proc. R. Soc. A 469: 20120348. (doi:10.1098/rspa.2012.0348).
  * [30] Strogatz SH. 1994 _Nonlinear Dynamics and Chaos_. Reading, US: Addison-Wesley.
  * [31] Thomas R. 1980 On the Relation Between the Logical Structure of Systems and Their Ability to Generate Multiple Steady States or Sustained Oscillations in Numerical Methods in the Study of Critical Phenomena. 180–193. Berlin, Germany: Springer (doi:10.1007/978-3-642-81703-8 24).
  * [32] Dawes JHP et al. $D_{n}$-symmetric dynamics under weak symmetry breaking. In preparation.
  * [33] Lloret-Fernández C, Maicas M, Mora-Martínez C, Artacho A, Jimeno-Martín Á, Chirivella L, Weinberg P, Flames N. 2018 A transcription factor collective defines the HSN serotonergic neuron regulatory landscape. eLife 7, e32785. (doi: 10.7554/eLife.32785).
  * [34] Jutras-Dubé L, El-Sherif E, François P. 2020 Geometric models for robust encoding of dynamical information into embryonic patterns. eLife 9, e55778. (doi:10.7554/eLife.55778).


# A Silicon Nitride Microring Modulator for High-Performance Photonic
Integrated Circuits

Venkata Sai Praneeth Karempudi  Ishan G Thakkar  Jeffrey Todd Hastings

###### Abstract

The use of the Silicon-on-Insulator (SOI) platform has been prominent for
realizing CMOS-compatible, high-performance photonic integrated circuits
(PICs). But in recent years, the silicon-nitride-on-silicon-dioxide (SiN-on-
SiO2) platform has garnered increasing interest as an alternative, because of
its several beneficial properties over the SOI platform, such as low optical
losses, high thermo-optic stability, broader wavelength transparency range,
and high tolerance to fabrication-process variations. However, SiN-on-SiO2
based active devices, such as modulators, are scarce and lack in desired
performance due to the absence of free-carrier-based activity in the SiN
material and the complexity of integrating other active materials with SiN-on-
SiO2 platform. This shortcoming hinders the SiN-on-SiO2 platform for realizing
active PICs. To address this shortcoming, in this article, we demonstrate a
SiN-on-SiO2 microring resonator (MRR) based active modulator. Our designed MRR
modulator employs an Indium-Tin-Oxide (ITO)-SiO2-ITO thin-film stack as the
active upper cladding and leverages the free-carrier assisted, high-amplitude
refractive index change in the ITO films to affect a large electro-refractive
optical modulation in the device. Based on the electrostatic, transient, and
finite difference time domain (FDTD) simulations, conducted using photonics
foundry-validated tools, we show that our modulator achieves 450 pm/V
resonance modulation efficiency, $\sim$46.2 GHz 3-dB modulation bandwidth, 18
nm free-spectral range (FSR), 0.24 dB insertion loss, and 8.2 dB extinction
ratio for optical on-off-keying (OOK) modulation at 30 Gb/s.

††journal: ol

##  1 Introduction

Driven by the rise of CMOS-compatible processes for fabricating photonic
devices, photonic integrated circuits (PICs) are inexorably moving from the
domain of long-distance communications to chip-chip and even on-chip
applications. It is common for the PICs to incorporate optical modulators to
enable efficient manipulation of optical signals, which is a necessity for the
operation of active PICs. Recent advances in the CMOS-compatible silicon-on-
insulator (SOI) photonic platform has fundamentally improved the applicability
of SOI PICs [1], [2], [3]. But in the last few years, the silicon-nitride-on-
silicon-dioxide (SiN-on-SiO2) platform has gained tremendous attention for
realizing PICs. This is because the SiN-on-SiO2 platform has several
advantages over the SOI platform. Compared to silicon (Si), the SiN material
has a much broader wavelength transparency range (500nm-3700nm), smaller
thermo-optic coefficient, and lower refractive index [4]. The lower refractive
index of SiN means that SiN offers smaller index contrast with SiO2 compared
to Si. This in turn makes the SiN-on-SiO2 based monomode passive devices
(e.g., waveguides, microring resonators (MRRs)) less susceptible to (i)
propagation losses due to the decreased sensitivity to edge roughness [5], and
(ii) aberrations in the realized device dimensions due to fabrication-process
variations [4]. In addition, the smaller thermo-optic coefficient of SiN makes
it possible to design nearly athermal photonic devices using SiN [6].
Moreover, SiN devices and circuits exhibit increased efficiency of nonlinear
parametric processes compared to Si [7].

Despite these favorable properties of the SiN-on-SiO2 platform, SiN-on-SiO2
based active devices such as modulators are scarce and lack in free spectral
range (FSR), modulation bandwidth, and modulation efficiency [8]. The lack in
efficiency is because of the lack of the free-carriers based activity in the
SiN material and the general difficulty of incorporating other active
materials with the SiN-on-SiO2 platform. This in turn limits the use of the
SiN-on-SiO2 platform to realizing only passive PICs. To overcome this
shortcoming, there is impetus to heterogeneously integrate active photonic
materials and devices with SiN-on-SiO2 passive devices (e.g., [9, 10, 11, 12,
13, 14, 15, 16, 17]). When such efforts of integrating electro-optically
active materials with the SiN-on-SiO2 platform come to fruition, it will be
possible to design high-performance and energy-efficient SiN-on-SiO2 based
active and passive PICs.

Different from such prior efforts, in this article, we demonstrate the use of
the high-amplitude electro-refractive activity of Indium-Tin-Oxide (ITO) thin
films to realize a SiN-on-SiO2 based optical on-off-keying (OOK) modulator. We
show, based on the electrostatic, transient, and finite difference time domain
(FDTD) simulations conducted using the photonics foundry-validated tools from
Lumerical/Ansys, that our modulator achieves 450 pm/V resonance modulation
efficiency, $\sim$46.2 GHz 3-dB modulation bandwidth, 18 nm free-spectral
range (FSR), 0.24 dB insertion loss, and 8.2 dB extinction ratio for optical
OOK modulation at 30 Gb/s. Based on the obtained simulation results, we
advocate that our modulator can achieve better performance compared to the
existing SiN modulators from prior works.

##  2 Related Work and Motivation

A plethora of Si and SiN based integrated optical modulator designs have been
formulated in prior works [18]. But among these modulator designs, MRR based
modulators have gained widespread attention due to their high wavelength
selectivity, compact size, and compatibility for cascaded dense wavelength
division multiplexing (DWDM). Recently several MRR based SiN-on-SiO2
modulators have also been demonstrated (e.g., [10, 8, 11, 12, 13, 14, 15, 16,
17]). In [10], a graphene integrated electro-optic SiN MRR modulator has been
reported. In [12], a hybrid SiN-LiNbO3 platform based racetrack resonator
modulator has been presented. Similarly, SiN modulators based on lead
zinconate titanate and zinc oxide/zinc sulphide as active materials are
demonstrated in [8] and [13]. In [11], a SiN modulator that achieves tuning
via photo-elastic effect has been demonstrated. Compared to these modulator
designs from prior works, we present a different, ITO-based electro-refractive
SiN-on-SiO2 modulator that achieves substantially better modulation bandwidth,
efficiency, and FSR.

Figure 1: (a) Top view, (b) Cross-sectional view (along AA’) of our SiN-on-
SiO2 MRR modulator.

##  3 Design of our SiN-on-SiO2 Modulator

In this section, firstly we describe the structure and operating principle of
our modulator design. Then, we discuss the characterization results for our
modulator that we have obtained through photonics foundry-validated
simulations. We also compare our modulator with several SiN based MRR
modulators from prior works, in terms of modulation bandwidth, modulation
efficiency, and FSR.

Figure 2: Transmission spectra of our modulator.

###  3.1 Structure and Operating Principle

Fig. 1(a) and Fig. 1(b), respectively, show the top-view and cross-sectional
schematic of our SiN-on-SiO2 MRR modulator. The active region in the upper
cladding of the modulator consists of a stack of two indium tin oxide (ITO)
thin films with a silicon dioxide (SiO2) thin film in between (an ITO-SiO2-ITO
stack). From Fig. 1(b), we have a 300 nm thick SiN-based MRR waveguide, two 10
nm thick ITO films, and 15 nm thick SiO2 layer. Upon applying voltage across
the ITO-SiO2-ITO stack (through the Au pads shown in Fig. 1(a)), free carriers
accumulate in the ITO films at the ITO-SiO2 interfaces for up to 5 nm depth in
the ITO films [1], making these accumulation regions in the ITO films high-
carrier-density active regions. In these regions, a free-carriers-assisted,
large-amplitude modulation in the permittivity and refractive index of the ITO
material has been previously reported [1]. We evaluate this free-carriers
based index modulation in the ITO films using the Drude-Lorentz model from
[19]. Accordingly, as the carrier concentration in the ITO accumulation
regions increases, the refractive index of the ITO films decreases. As a
result, the effective refractive index of our SiN-on-SiO2 modulator design
from Fig. 1 also decreases, causing a blue shift in its resonance wavelength
that in turn causes a transmission modulation at the through port of the
modulator. The electro-refractive activity of our SiN-on-SiO2 MRR modulator
design is confined only in the ITO-SiO2-ITO cladding. This is different from
the Si-SiO2-ITO capacitor based MRR modulator from [20], which has the
electro-refractive as well as electro-absorptive activities in both its Si-
based MRR core and SiO2-ITO based cladding.

Table 1: Free-carrier concentration (N), real index (Re($\eta_{ITO}$)), and imaginary index (Im($\eta_{ITO}$)) for the ITO accumulation layer in our modulator. The real and imaginary effective index (Re($\eta_{eff}$), Im($\eta_{eff}$)), operating voltage (V), and induced resonance shift ($\Delta\lambda_{r}$) for our modulator. | N
---
($cm^{-3}$)
| Re
---
($\eta_{ITO}$)
| Im
---
($\eta_{ITO}$)
| Re
---
($\eta_{eff}$)
| Im
---
($\eta_{eff}$)
V |  | $\Delta\lambda_{r}$
---
(pm)
1 × $10^{19}$ | 1.9556 | 0.0100 | 1.9735 | 0.0001 | 0 | 0
5 × $10^{19}$ | 1.9111 | 0.0403 | 1.9724 | 0.0003 | 1.8 | 830
9 × $10^{19}$ | 1.8667 | 0.0896 | 1.9712 | 0.0006 | 3.7 | 1580
13 × $10^{19}$ | 1.8222 | 0.1289 | 1.9701 | 0.0011 | 5.5 | 2470
17 × $10^{19}$ | 1.7778 | 0.1582 | 1.9692 | 0.0017 | 7.3 | 3210
20 × $10^{19}$ | 1.7333 | 0.1874 | 1.9680 | 0.0022 | 9.2 | 4000

###  3.2 Simulations Based Characterization

We performed electrostatic simulations of our ITO-SiO2-ITO thin-film stack
based SiN-on-SiO2 modulator in the CHARGE tool of DEVICE suite from Lumerical
[21], to evaluate the required voltage levels across the Au pads (Fig. 1(a))
for achieving various free-carrier concentrations in the ITO films. Then,
based on the Drude-Lorentz dispersion model from [19], we extracted the
corresponding ITO index change values for various free-carrier concentrations.
These results are listed in Table 1. Using these index values from Table 1, we
modeled our MRR modulator in the MODE tool from Lumerical [21] for finite-
difference-time-domain (FDTD) and finite-difference eigenmode (FDE) analysis.
For this analysis, we used the Kischkat model [22] of stoichiometric silicon
nitride to model the MRR device. From this analysis, we extracted the
effective index change and transmission spectra of our modulator (shown in
Table 1 and Fig. 2 respectively) at various applied voltages for the operation
around 1.6 $\mu$m wavelength (L-band). From Fig. 2, our modulator achieves up
to 4 nm resonance shift upon applying 9.2 V across the thin-film stack, which
renders the resonance tuning (modulation) efficiency of $\sim$450 pm/V. This
is a crucial outcome as our MRR modulator has relatively very low overlap
between the optical mode and free-carrier perturbation (only 10% of the guided
optical mode overlaps with the ITO-based upper cladding) compared to the
silicon ITO-based modulators (e.g., [20]).

Figure 3: Cross-sectional electric-field profiles of the fundamental TE mode
evaluated at the coupling section (along BB’in Fig. 1(a)) ((a)-(c)), across
the rim (along AA’ in Fig. 1) ((d)-(f)), and at the through port of our SiN-
on-SiO2 MRR modulator ((g)-(i)), for three different free-carrier
concentrations of ITO (Table 1) namely 1$\times$1019 cm-3 (for (a),(d),(g)),
9$\times$1019 cm-3 (for (b),(e),(f)), and 17$\times$1019 cm-3 (for
(c),(f),(i)), using the variational FDTD (varFDTD) solver [21].

Fig. 3 illustrates the cross-sectional electric-field profiles of the
fundamental TE mode evaluated for three different free-carrier concentrations
of ITO, namely 1$\times$1019 cm-3, 9$\times$1019 cm-3, and 17$\times$1019
cm-3, at three different cross-sectional regions of our MRR modulator. To
evaluate these profiles, we used the variational FDTD (varFDTD) solver of the
MODE tool of the DEVICE suite from ANSYS/Lumerical. In fact, Fig. 3 shows a
3$\times$3 grid of field profiles. Each row in this grid corresponds to field
profiles collected for a particular cross-sectional region of our MRR
modulator across different free-carrier concentrations of ITO. Similarly, each
column in the grid corresponds to field profiles collected for a particular
free-carrier concentration of ITO across three different regions of the
modulator.

As per the discussion in the previous subsection, the increase in the free-
carrier concentration in the ITO layers caused due to the increase in the
applied bias across the ITO-SiO2-ITO stack, decreases the effective index of
our modulator. This in turn induces a blue shift in the resonance wavelength
of our modulator. Due to this blue shift in the resonance wavelength, the
amount of optical power coupled from the input port into the MRR cavity at the
coupling region decreases. This can be clearly observed from the field
profiles collected at the coupling region along BB’; as the free-carrier
concentration increases from Fig. 3(a) to Fig. 3(c), the intensity of the
coupled field in the MRR at the cross-section BB’ also decreases. The decrease
in the coupled field intensity at BB’ naturally results in the decrease of the
steady-state field intensity inside the MRR waveguide. As a result, at the
cross-section AA’, the field intensity can be observed to decrease with the
increase in the free-carrier concentration in the ITO layers, as we move from
Fig. 3(d) to Fig. 3(f) in the middle row of Fig. 3. Atop the steady-state
field intensity inside the MRR cavity, the field intensity at the through port
(hence, the output optical power at the through port) of the MRR also
decreases naturally with the increase in the free-carrier concentration. This
can be observed in the bottom row of Fig. 3. The modulation of the optical
output power at the through port with the change in the free-carrier
concentration in the ITO layers corroborates the electro-refractive activity
of our modulator.

In addition, as we move from the top row (coupling region field profiles) to
the bottom row (through port field profile) within each column in Fig. 3, the
field intensity slightly decreases. This provides evidence that, for each
column (i.e., for each specific free-carrier concentration), the optical field
intensity undergoes optical-loss-induced attenuation as the light waves travel
along the propagation path from the coupling region (top row) to the through
port (bottom row).

Figure 4: Optical eye diagrams for (a) 30 Gb/s and (b) 55 Gb/s OOK inputs to
our modulator.

Further, from the spectra in Fig. 2, we evaluate the FSR of our modulator to
be $\sim$18 nm. We evaluated (using the Lumerical MODE tool) the insertion
loss and loaded Q-factor of our modulator to be $\sim$0.235 dB and $\sim$2000
respectively. We also evaluated the capacitance density of the ITO thin-film
stack covering the MRR rim (using the Lumerical CHARGE tool) to be $\sim$2.3
fF$/\mu$m2 for the 15 nm thick SiO2 layer. Moreover, we modeled our modulator
in Lumerical INTERCONNECT, to simulate optical eye diagrams for the modulator
at 30 Gb/s and 55 Gb/s operating bitrates (Fig. 4). As evident (Fig. 4), our
modulator can achieve 8.2 dB extinction ratio for OOK modulation at 30 Gb/s
bitrate.

Table 2: Modulation bandwidth (MB) (optical (O) and Electrical (E)), modulation efficiency (ME), FSR and energy efficiency (EE) corresponding to various SiN based MRR modulators (modulator type (MT)) from prior works obtained from simulations (*) and experiments, compared with our simulated SiN-on-SiO2 MRR modulator. MT | MB |  ME (pm/V) |  FSR (nm) |  EE (pJ/bit)
---|---|---|---|---
| O-MB
---
(GHz)
| E-MB
---
(GHz)
[14] | 0.06 | 0.02 | 1.6 | 4×10-3 | 1×10-3
[15] | 1.7 | N/A | 0.01* | 0.58 | N/A
[16] | 0.01 | 7.9 | 1 | 113 | 3.7×104
[17] | 0.03 | 0.03 | 1.6 | 0.3 | 1×10-3
[10] | 161 | 30 | 100* | 4.7 | 0.8
[8] | 87 | 35.6 | 67* | 1.74 | N/A
[11] | 0.19 | 1.3×10-3 | 137.5* | 2.2* | 0.11
[12] | 25 | N/A | 2.9 | 0.21* | N/A
[13] | 1.55 | 5.9 | 0.2 | 0.6 | 53
[23] | 1.77 | N/A | 5.8 | 0.4* | N/A
| Ours
---
93.62* | 53.1* | 450* | 18* | 1.4*

###  3.3 Comparison and Discussion

Table 2 shows a comparison of our SiN-on-SiO2 MRR modulator with the
simulation (marked as *) and fabrication based best-performing nine SiN MRR
modulators from prior works ([10]-[13],[14]-[17],[23]), in terms of five key
attributes namely optical modulation bandwidth (O-MB), electrical modulation
bandwidth (E-MB), modulation efficiency (ME), FSR, and energy-efficiency (EE).
The SiN MRR modulator in [10] achieves higher O-MB compared to the other SiN
MRR modulators (Table 2) and our modulator. In contrast, our modulator
achieves higher E-MB compared to the other SiN MRR modulators (Table 2). We
have also evaluated that our modulator achieves the best effective MB of
$\sim$46.2 GHz compared to all other SiN MRR modulators, based on the formula
of effective MB from [24]. Due to its superior effective MB of $\sim$46.2 GHz,
our modulator can be easily operated at $>$15 Gb/s bitrate to enable ultra-
high-speed (potentially beyond Tb/s) DWDM-based PICs while ensuring minimal
power-penalty from crosstalk [25]. Moreover, our modulator achieves higher ME
compared to other SiN MRR modulators (Table 2). However, in terms of FSR, SiN
MRR modulator demonstrated in [16] achieves higher FSR compared to the other
SiN MRR modulators in Table 2 including our modulator. Nevertheless, our
modulator consumes the energy of 1.4 pJ/bit which is significantly better than
the energy consumption of the modulator from [16]. Its high energy efficiency
and competent FSR of 18 nm make our modulator a favorable candidate for
designing high-bandwidth and energy-efficient DWDM-based photonic
interconnects for datacenter-scale as well as chip-scale computing and
communication architectures.

Further, although ITO is not available in the CMOS process flow, it can be
deposited at relatively low temperatures (less than 300°C) on top of the back-
end-of-line (BEOL) metal layers of CMOS chips, in an independent manner
without interfering with or contaminating the CMOS front-end-of-line (FEOL)
and BEOL processes. This makes our SiN-on-SiO2 modulator an excellent choice
for implementing optical interconnect PICs on silicon interposers, to enable
ultra-high-bandwidth inter-chiplet communication in emerging multi-chiplet
systems [26].

##  4 Conclusion

We have demonstrated an ITO-based SiN-on-SiO2 MRR modulator, which consists of
a stack of ITO-SiO2-ITO thin films as the active upper cladding of the SiN MRR
core. This active upper cladding of our modulator leverages the free-carrier
assisted, high-amplitude refractive index change in the ITO films to affect a
large electro-refractive optical modulation in the device. To evaluate the
performance of our SiN-on-SiO2 MRR modulator, we performed electrostatic,
transient, and FDTD simulations using the foundry-validated Ansys/Lumerical
tools. Based on these simulations, our modulator achieves superior performance
with 450 pm/V modulation efficiency, $\sim$46.2 GHz 3-dB modulation bandwidth,
18nm FSR, 0.24 dB insertion loss, and 8.2 dB extinction ratio for OOK
modulation at 30 Gb/s. This excellent performance of our SiN-on-SiO2 MRR
modulator demonstrates its potential to enhance the performance and energy
efficiency of SiN-on-SiO2 based PICs of the future.

Disclosures The authors declare no conflicts of interest

## References

  * [1] L. Chrostowski, J. Flueckiger, C. Lin, M. Hochberg, J. Pond, J. Klein, J. Ferguson, and C. Cone, SPIE SPOIC XVI 8989, 83 (2014).
  * [2] W. Bogaerts, X. Chen, M. Wang, I. Zand, H. Deng, L. Van Iseghem, A. Ribeiro, A. D. Tormo, and U. Khan, IEEE IPC pp. 1–2 (2020).
  * [3] N. C. Harris, J. Carolan, D. Bunandar, M. Prabhu, M. Hochberg, T. Baehr-Jones, M. L. Fanto, A. M. Smith, C. C. Tison, P. M. Alsing, and D. Englund, Optica 5, 1623 (2018).
  * [4] Q. Wilmart, H. El Dirani, N. Tyler, D. Fowler, S. Malhouitre, S. Garcia, M. Casale, S. Kerdiles, K. Hassan, and C. Monat, Applied Sciences 9, 255 (2019).
  * [5] J. F. Bauters, M. J. Heck, D. John, D. Dai, M.-C. Tien, J. S. Barton, A. Leinse, R. G. Heideman, D. J. Blumenthal, and J. E. Bowers, Optica OE 19, 3163 (2011).
  * [6] G. Gao, D. Chen, S. Tao, Y. Zhang, S. Zhu, X. Xiao, and J. Xia, Optica OE 25, 12260 (2017).
  * [7] J. Levy, Cornell Theses and Dissertations (2011).
  * [8] K. Alexander, J. P. George, J. Verbist, K. Neyts, B. Kuyken, D. Van Thourhout, and J. Beeckman, Nature communications 9, 1 (2018).
  * [9] J. Goyvaerts, A. Grabowski, J. Gustavsson, S. Kumari, A. Stassen, R. Baets, A. Larsson, and G. Roelkens, Optica 8, 1573 (2021).
  * [10] C. T. Phare, Y.-H. Daniel Lee, J. Cardenas, and M. Lipson, Nature Photonics 9, 511 (2015).
  * [11] W. Jin, R. G. Polcawich, P. A. Morton, and J. E. Bowers, Optica OE 26, 3174 (2018).
  * [12] A. N. R. Ahmed, S. Shi, A. J. Mercante, and D. W. Prather, Optica OE 27, 30741 (2019).
  * [13] A. Hermans, M. Van Daele, J. Dendooven, S. Clemmen, C. Detavernier, and R. Baets, Optica OL 44, 1112 (2019).
  * [14] J. Wang, K. Liu, M. W. Harrington, R. Q. Rudy, and D. J. Blumenthal, Optica OE 30, 31816 (2022).
  * [15] A. N. R. Ahmed, S. Shi, M. Zablocki, P. Yao, and D. W. Prather, Optica OL 44, 618 (2019).
  * [16] J. Liu, H. Tian, E. Lucas, A. S. Raja, G. Lihachev, R. N. Wang, J. He, T. Liu, M. H. Anderson, and W. Weng, Nature 583, 385 (2020).
  * [17] J. Wang, K. Liu, M. W. Harrington, R. Q. Rudy, and D. J. Blumenthal, IEEE OFC pp. 1–3 (2022).
  * [18] A. Rahim, A. Hermans, B. Wohlfeil, D. Petousi, and B. Kuyken, SPIE Advanced Photonics 3, 024003 (2021).
  * [19] Z. Ma, Z. Li, K. Liu, C. Ye, and V. J. Sorger, Nanophotonics 4, 198 (2015).
  * [20] E. Li, B. A. Nia, B. Zhou, and A. X. Wang, IEEE OI pp. 1–2 (2019).
  * [21] Lumerical, “https://www.ansys.com/products/photonics,” .
  * [22] J. Kischkat, S. Peters, B. Gruska, M. Semtsiv, M. Chashnikova, M. Klinkmüller, O. Fedosenko, S. Machulik, A. Aleksandrova, and G. Monastyrskyi, Optica AO 51, 6789 (2012).
  * [23] Z. Zhao, Z. Zhang, J. Li, Z. Shang, G. Wang, J. Yin, H. Chen, K. Guo, and P. Yan, Optica OL 47, 949 (2022).
  * [24] O. Dubray, A. Abraham, K. Hassan, S. Olivier, D. Marris-Morini, L. Vivien, I. O’Connor, and S. Menezo, IEEE JSTQE 22, 89 (2016).
  * [25] M. Bahadori, S. Rumley, H. Jayatilleka, K. Murray, N. A. Jaeger, L. Chrostowski, S. Shekhar, and K. Bergman, IEEE JLT 34, 4043 (2016).
  * [26] Y. Thonnart, S. Bernabe, J. Charbonnier, C. Bernard, D. Coriat, C. Fuguet, P. Tissier, B. Charbonnier, S. Malhouitre, D. Saint-Patrice, M. Assous, A. Narayan, A. Coskun, D. Dutoit, and P. Vivet, IEEE/ACM DATE pp. 1456–1461 (2020).

ref