1. Introduction

Transcription factors play a vital role in controlling the levels of proteins and mRNAs within cells, and are involved in many key processes such as cell-cycle regulation, inflammation, meiosis, apoptosis and the heat shock response [20]. Such systems are often referred to as gene regulatory networks (GRNs). Those transcription factors which down-regulate (repress/suppress) the rate of gene transcription do so via negative feedback loops, and such intracellular negative feedback systems are known to exhibit oscillations in protein and mRNA levels [18]. It is now known that GRNs contain a small set of recurring regulation patterns, commonly referred to as network motifs [25], which can be thought of as recurring circuits of interactions from which complex GRNs are built. When studied from a theoretical point of view, such intracellular systems of interacting network motifs have been referred to as synthetic genetic oscillators or synthetic gene oscillators [9,29,30].

The Hes1 system is a simple example of a GRN which possesses a single negative feedback loop and benefits from having been the subject of numerous biological experiments [14,16,17,18,19] as well as attention from a more theoretical perspective via mathematical models. Mathematical modelling of Hes1 can be traced back 50 years to the seminal paper of Goodwin [11], followed shortly after by the paper of Griffith [13]. The introduction of time delays by Tiana et al. [39] led to a period of activity of the study of delay-differential equation models of GRNs, including the Hes1 system, the p53-Mdm2 system and the NF-$\kappa$B system [39,15,27,26].

Early spatial models of theoretical intracellular systems were pioneered in the 1970s by Glass and co-workers [10,33] and again in the 1980s by Mahaffy and co-workers [1,23,24], where the focus was on generic systems with one-dimensional models. This approach has been extended by Naqib et al. [28]. More recently, a two-dimensional spatial model of molecular transport inside a cell was formulated by Cangiani and Natalini [2] and this general approach was adopted by Sturrock et al. [36] to formulate and study a spatio-temporal model of the Hes1 GRN considering diffusion of the protein and mRNA. This model was then later extended to account for transport across the nuclear membrane and directed transport via microtubules [37] and finally to develop a spatial stochastic model of the Hes1 system [35]. Other papers adopting an explicitly spatial approach include those of Szyma\'{n}ska et al. [38] focussing on the role of transport via the microtubules and Clairambault and co-workers [5,6,7,8] focussing on the p53 system. The need to investigate the role that spatial components play has been bolstered by the results of experimentalists which show the importance of the spatial location of genes (cf. [21]).

In this paper we first formulate and study a mathematical model of the Hes1 transcription factor -a canonical GRN consisting of a single negative feedback loop between the Hes1 protein and its mRNA. A key feature of the Hes1 GRN (and other negative feedback systems) is the existence of oscillatory solutions in the levels of mRNA and protein. Following on from the results of [3], where it was proven that the diffusion coefficients of the mRNA and protein were Hopf bifurcation parameters, we explore the effect of varying the distance between gene sites (transcription) and protein production sites (translation) on the oscillations. This approach is then generalized to explore the effect of varying the distances between gene and protein production sites on the spatio-temporal dynamics of synthetic GRNs, such as repressilators and activator-repressor oscillators.

---

2. The Hes1 system: The canonical gene regulatory network

The Hes1 protein is a member of the family of basic helix-loop-helix (bHLH) transcription factors. Hes1 is known to repress the transcription of its own gene through direct binding to regulatory sequences in the Hes1 promoter [14], and hence may be termed the canonical gene regulatory network. The Hes1 GRN is involved in key processes during embryogenesis [16] and plays a central role in the timing of somitogenesis [14]. It also plays a role in the development of certain human cancers (where it can become deregulated [31]) and in the differentiation of embryonic stem cells [17,18,19].

In this section we present a spatial model of the Hes1 GRN based on a previous model originally formulated by Sturrock et al. [36,37]. The model consists of a system of coupled nonlinear partial differential equations describing the temporal and spatial dynamics of the concentration of hes1 mRNA, $m(x, t)$, and Hes1 protein, $p(x, t)$, in a 1-dimensional domain. and accounts for the processes of transcription (mRNA production) and translation (protein production). Transcription is assumed to occur in a small region of the domain representing the gene site. Both mRNA and protein also diffuse and undergo linear decay. The (non-dimensionalised) model is given as:

$\frac{\partial m }{\partial t} = D_m \frac{\partial^2 m }{\partial x^2} + \dfrac{\alpha_m}{1+p^h}\delta^\varepsilon _{x_m}(x) - \mu_m \, m \; \text{ in }\; (0, T)\times (-1, 1), \\ \frac{\partial p}{\partial t} = D_p \frac{\partial^2 p }{\partial x^2}+ \alpha_p \, m \, \delta^\varepsilon _{x_p}(x) - \mu_p \, p \; \text{ in }\; (0, T)\times (-1, 1), \\ \frac{\partial m(-1, t)}{\partial x} =\frac{\partial m(1, t)}{\partial x} =0, \frac{\partial p(-1, t)}{\partial x} = \frac{\partial p(1, t)}{\partial x} =0 \; \text{ in }\; (0, T), \\ m(0, x)= m_0(x), p(0, x) = p_0(x) \; \text{ in }\; (-1, 1)\;,$

(1)

where $D_m, D_p$, $\alpha_m$, $\alpha_p$, $\mu_m$ and $\mu_p$ are positive constants (the diffusion coefficient, transcription rate, translation rate and decay rates of hes1 mRNA and Hes1 protein respectively). The nonlinear reaction term ${1}/{(1+p^h)}$, with $h \geq 2$, is a Hill function, modelling the suppression of mRNA production by the protein (negative feedback). The point $x_m$ is the position of the centre of the gene site (where transcription occurs) and similarly the point $x_p$ is the position of the centre of protein production (where translation occurs). We use a Dirac approximation of the $\delta$-distribution function located at the gene site and protein production sites $x_i$, where $i=\{m, \pm p\}$, such that

$\delta_{x_i}^{\varepsilon}(x) = \left\{ \begin{array}{cc} \dfrac{1}{2\varepsilon}(1+\cos(\pi(x-x_i)/\varepsilon)), & |x-x_{i}| < \varepsilon, \\ & \\ 0, & |x-x_{i}| \ge \varepsilon, \end{array}\right.$

(2)

with $\varepsilon > 0$ a small parameter. In fact, unless otherwise stated, for the remainder of the paper we take $x_m = 0$. We consider a symmetric 1D interval [-1, 1] which may be thought of as the cross section of a cell, with the mRNA being produced in the nucleus and the protein being produced in the cytoplasm. The localisation of protein production around a single point $x_p$ is a simplification since translation occurs via the ribosomes which are located throughout the cytoplasm. However, this provides us with a useful parameter (i.e. the distance between $x_m$ and $x_p$) when assessing whether or not oscillations occur, and in a later section we do investigate the effect of a uniform protein production in line with the original models of Sturrock et al. [36,37].

Unless otherwise stated, the parameter values we will use throughout the paper are as follows:

$D_m = D_p = 0.00075; \alpha_m = 1.0; \alpha_p = 2.0; \mu = 0.03; h = 5; \varepsilon = 0.01,$

which we refer to as the baseline parameter set $\mathbb{P}$. Full details of the scalings and non-dimensionalisation can be found in [36,37].

---

2.1. Computational simulation results

Numerical simulations of system (1) in [36] demonstrated the existence of oscillatory solutions as observed experimentally, with the indication that the periodic orbits arose from supercritical Hopf bifurcations at two critical values of the diffusion coefficients $D_m, D_p$. Recent work by Chaplain, Ptashnyk and Sturrock [3] analysing this model proved rigorously that the diffusion coefficients of the protein/mRNA act as bifurcation parameters, thus showing that the spatial movement of the molecules alone is sufficient to cause the oscillations. In this section we present computational simulation results exploring the effect of varying the distance between the gene site (where transcription occurs) and the protein production sites (where translation occurs) on the spatio-temporal dynamics of the system.

Figure 1 shows the computational results of a numerical simulation of (1) where the hes1 mRNA gene-site (where transcription occurs) is located at $x_{m}=0.0$ and where the Hes1 protein production sites (where translation occurs) are located at $x_{p}=\pm0.5$. The plots show sustained oscillations of both hes1 mRNA and Hes1 protein in space and time. For this and all other computational results we use the pdepe solver in Matlab which uses a method of lines approach with a stiff ODE solver (cf. [34]).

Figure 2 shows the computational results of a numerical simulation of (1) where the locations of the Hes1 protein production sites are varied relative to the fixed hes1 mRNA gene-site at $x=0$. The plots show the total concentrations of hes1 mRNA and Hes1 protein over time and demonstrate that if the Hes1 protein production sites are either too close or too far away from the hes1 mRNA gene-site, then oscillations are lost.

---

2.1.1. Modification of the protein production term

As mentioned in Section 2, translation (i.e. protein production) occurs via the ribosomes which are located throughout the cytoplasm. To model this more accurately and to make comparisons with that of [36,37] we now consider two scenarios where the protein production term $\alpha_p \, m \, \delta^\varepsilon _{x_p}(x)$ is modified to take note of this fact. In the first case we modify the Dirac approximation by increasing the parameter $\varepsilon$ (now taken to be $\varepsilon =0.25$), such that translation occurs over a wider region in the cytoplasm. In the second case we replace the Dirac approximation with a Heaviside function which means that translation occurs uniformly in a region $|x|\ge x_p$, i.e. the equation for $p$ in (1) becomes

$\frac{\partial p}{\partial t} = D_p \frac{\partial^2 p }{\partial x^2}+ \alpha_p \, m \, H_{x_p} (x) - \mu_p \, p \;\;\; \text{ in }\; (0, T)\times (-1, 1),$

(3)

where

$H_{x_p} (x) = \left\{ \begin{array}{cl} 1, & | x | > x_p \\ 0, & | x | < x_p . \end{array}\right.$

(4)

The results for these simulations are shown in Figures 3 and 4 respectively. We note that they are qualitatively the same as those shown in Figure 2 and that our conclusion, that protein production must be neither too close to nor too far apart from mRNA production, remains valid.

---

3. Synthetic GRNs - n-gene repressilators

In this section we consider a GRN known as a repressilator which can be defined as a system of interactions between multiple genes/proteins in which the protein of any given gene inhibits the production of the mRNA for the subsequent gene [9,29,30]. In fact, the Hes1 system discussed in the previous section can be viewed as a one-gene repressilator and so we base our synthetic multi-gene repressilator on the Hes1 system structure. As such the equations in 1D are taken to be:

$\begin{split} \dfrac{\partial m_i}{\partial t} = D_{m_i} \dfrac{\partial^2 m_i}{\partial x^2}+\dfrac{\alpha_m}{1+p_j^h}\delta_{x_{m_i}}^{\varepsilon}(x)-\mu_mm_i, \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\dfrac{\partial p_i}{\partial t} = D_{p_i} \dfrac{\partial^2 p_i}{\partial x^2}+\alpha_p m_i \delta_{x_{p_i}}^{\epsilon}(x)-\mu_pp_i, \end{split}$

(5)

where $i={1, 2, 3, \ldots n}$ and, since the repression of mRNA comes from the preceding protein in the system, $j=\{n, 1, 2, 3, \ldots (n-1)\}$ for an $n$-gene system. As for the Hes1 system we use a Dirac approximation of the $\delta$-distribution function, this time located at the gene-sites $x_{m_i}$ (transcription) and protein production site $x_{p_i}$ (translation), with $i$ as above. The boundary conditions and initial conditions are, as before, such that:

$\begin{split} \dfrac{\partial m_i(-1, t)}{\partial x} = \dfrac{\partial m_i(1, t)}{\partial x} =0, \\ \dfrac{\partial p_i(-1, t)}{\partial x} = \dfrac{\partial p_i(1, t)}{\partial x} =0, \end{split}$

(6)

$m_i(x, 0)=0,\;\;\;\;\;\; p_i(x, 0)=0.$

(7)

---

3.1. Computational simulation results

In this section we present computational simulation results from a 3-gene repressilator system i.e. system (5) with $i=1, 2, 3$ and $j=3, 1, 2$. We explore the effect of varying the distance between the gene site (where transcription occurs) and the protein production sites (where translation occurs) on the spatio-temporal dynamics of the system.

Figure 5 shows the computational results of a numerical simulation where the three mRNA gene-sites are located at $x_{m1}=0.0$, $x_{m2}=\pm0.2$, $x_{m3}=\pm0.4$, and the three protein production sites are located at $x_{p1}=\pm0.5$, $x_{p2}=\pm0.7$ and $x_{p3}=\pm0.9$. The plots show sustained oscillations of all three species mRNA and protein in space and time.

Figure 6 shows the computational results of a numerical simulation of (1) where the locations of the gene and protein production sites are varied relative to each other. The plots show the total concentrations of mRNA and protein for species 1 over time and demonstrate that if distance between gene production sites and protein production sites is too great, then oscillations are lost (similar results for species 2 and 3).

Figure 7 shows the computational results of a numerical simulation of (1) where the locations of the gene sites are varied but protein production is permitted throughout the cytoplasm which is fixed at $|x|>0.5$ using a Heaviside function as discussed in Section 2.1.1 (see Equations (3) and (4) with $x_p=0.5$). The plots show the total concentrations of mRNA and protein for species 1 over time and demonstrate that the precise spatial aspects control the nature of oscillations.

---

4. Activator-repressor oscillators

The second type of synthetic system we consider is one in which there is both positive and negative feedback, such that one protein activates a second protein through positive feedback at the mRNA level (i.e. up-regulation of mRNA), while the second protein represses the first (as per previous section). We refer to this as an activator-repressor system. We first consider the simplest such system i.e. two species each with its own mRNA and protein, and we modify the system of equations given for the two-gene repressilator system by modifying the production term for the second species so that now its mRNA production is increased by the first protein. The equations in 1D are:

$\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\dfrac{\partial m_1}{\partial t} = D_{m_1} \dfrac{\partial^2 m_1}{\partial x^2}+\dfrac{\alpha_m}{1+p_2^h}\delta_{x_{m_1}}^{\varepsilon}(x)-\mu_mm_1\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\dfrac{\partial p_1}{\partial t} = D_{p_1} \dfrac{\partial^2 p_1}{\partial x^2}+\alpha_p m_1 \delta_{x_{p_1}}^{\varepsilon}(x)-\mu_pp_1\\ \dfrac{\partial m_2}{\partial t} = D_{m_2} \dfrac{\partial^2 m_2}{\partial x^2}+ \left[\alpha_m + \dfrac{\beta_m p_1^h}{1+p_1^h} \right] \delta_{x_{m_2}}^{\varepsilon}(x)-\mu_mm_2\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\dfrac{\partial p_2}{\partial t} = D_{p_2} \dfrac{\partial^2 p_2}{\partial x^2}+\alpha_p m_2 \delta_{x_{p_2}}^{\varepsilon}(x)-\mu_pp_2$

(8)

where the $m_i(x, t)$ and $p_i(x, t)$ are, as before, the concentrations of mRNA and protein, respectively for genes $i=\{1, 2\}$. In the equation for $m_2$, the enhanced production term is given by:

$\alpha_m + \dfrac{\beta_m p_1^h}{1+p_1^h} \;,$

which is similar to that used by Sturrock et al. [36,37] in their model of the p53-Mdm2 system. The p53-Mdm2 system may be considered an activator-repressor system, but with the negative feedback being provided by Mdm2-enhanced p53 degradation via ubiquitination. In our system, the repression is provided directly by negative feedback to mRNA1, $m_1$, by protein2, $p_2$, as per the previous $n$-gene repressilator models.

The boundary conditions and the initial conditions are, as before,

$\begin{split} \dfrac{\partial m_i(-1, t)}{\partial x} = \dfrac{\partial m_i(1, t)}{\partial x} =0, \\ \dfrac{\partial p_i(-1, t)}{\partial x} = \dfrac{\partial p_i(1, t)}{\partial x} =0, \end{split}$

(9)

$m_i(x, 0)=0, \;\;\;\;\;\;\;\; p_i(x, 0)=0.s$

(10)

---

4.1. Computational simulation results

In this section we present some computational simulation results from the simplest activator-repressor system (8) with two species. Once again we explore the effect of varying the distance between the gene site (where transcription occurs) and the protein production sites (where translation occurs) on the spatio-temporal dynamics of the system. We take $\beta_m=10.0$ while all other parameters are as in baseline parameter set $\mathbb{P}$.

Figure 8 activator_repressor_spatial shows the computational results of a numerical simulation of (8) where the two mRNA gene-sites are located at $x_{m_i}=0.0, i=1, 2$ and where the two protein production sites are located at $x_{p_i}= \pm0.9, i=1, 2$. The plots show sustained oscillations of mRNA and protein for both species in space and time.

Figure 9 shows the computational results of a numerical simulation of (1) where the locations of the protein production sites are varied relative to the fixed mRNA gene-sites at $x=0$. The plots show the total concentrations of mRNA and protein for species 1 over time and demonstrate that if the protein production sites are too far away from the mRNA gene-site, then oscillations are lost (similar results for species 2).

---

5. Discussion

In this paper we have formulated and investigated spatio-temporal models of a range of gene regulatory networks, GRNs. The first model considered the canonical gene regulatory network, the Hes1 system, where the Hes1 protein is known to suppress hes1 mRNA production through a negative feedback. This system may also be viewed as a one-gene repressilator. This led us to then investigate the spatio-temporal dynamics of synthetic GRNs known as $n$-gene repressilators. These are cyclic feedback systems whereby each protein in the cycle suppresses the mRNA of the subsequent species, with the $n$-th protein suppressing the mRNA of the first species. Finally we examined the spatio-temporal dynamics of a two species activator-repressor system, where species one activated species two, which in turn repressed species one.

All such systems are known to have oscillatory dynamics in the concentration levels of the various species, but there have been relatively few spatio-temporal models (i.e. PDEs) of such systems. Inspired by the results of [3] where it was shown that the diffusion coefficient of the Hes1 system was a Hopf bifurcation parameter and caused the oscillations, and the work of [28] where similar results were obtained for an infinite domain, in this paper we have investigated the role of the location of the gene-site (where transcription occurs) and the protein production site (where translation occurs) on the oscillatory dynamics of synthetic GRNs. The results of our computational simulations have shown that the separation between mRNA and protein production sites for the simple Hes1 system must be optimised in order to achieve oscillations. If the distance between the sites is either too large or too small oscillations are lost. Similar results were obtained for the synthetic GRNs considered, both the $n$-gene repressilators and the activator-repressor system. Our computational simulation results, along with the results of [3,28] indicate the importance of spatial components when modelling GRNs. More generally, oscillatory spatio-temporal behaviour of molecules has also recently been identified as playing a role in the context of pattern formation (cf. [22]).

Future work in this area will consider a deeper study of more complex $n$-gene repressilators and activator-repressor systems, although preliminary results on such systems (e.g. $4-$, $5-$, $6-$gene repressilators) show qualitatively similar behaviour, as do models in 2-dimensional circular/elliptical and 3-dimensional spherical domains. Exploring these synthetic GRNs using spatial stochastic models (cf. [35]) is another avenue of investigation. The theoretical study of synthetic GRNs is very timely and makes certain predictions that can be used by experimentalists. Given the current level of interest in synthetic biology and the technological tools available to synthetic biologists [4,32,40], the findings in this paper indicate that experimentalists should take molecular movement into account when trying to design such systems.

---