Most probable dynamics of a genetic regulatory network under stable L\'evy noise

TL;DR

This paper investigates how non-Gaussian stable Le9vy noise influences the most probable trajectories of a genetic regulatory network, revealing parameter effects on state transitions and metastability in B. subtilis.

Contribution

It introduces a novel analysis of the MeKS network under Le9vy noise, highlighting how non-Gaussian noise parameters affect state transitions and metastable behaviors.

Findings

Optimal noise parameters shorten tipping time.

Different initial states evolve to metastable states.

Non-Gaussian noise influences transition pathways.

Abstract

Numerous studies have demonstrated the important role of noise in the dynamical behaviour of a complex system. The most probable trajectories of nonlinear systems under the influence of Gaussian noise have recently been studied already. However, there has been only a few works that examine how most probable trajectories in the two-dimensional system (MeKS network) are influenced under non-Gaussian stable L\'evy noise. Therefore, we discuss the most probable trajectories of a two-dimensional model depicting the competence behaviour in B. subtilis under the influence of stable L\'evy noise. On the basis of the Fokker-Planck equation, we describe the noise-induced most probable trajectories of the MeKS network from the low ComK protein concentration (vegetative state) to the high ComK protein concentration (competence state) under stable L\'evy noise. We demonstrate choices of the non-Gaussianity index $\alpha$ and the noise intensity $\epsilon$ which generate the ComK protein escape from the low concentration to the high concentration. We also reveal the optimal combination of both parameters $\alpha$ and $\epsilon$ making the tipping time shortest. Moreover, we find that different initial concentrations around the low ComK protein concentration evolve to a metastable state, and provide the optimal $\alpha$ and $\epsilon$ such that the distance between the deterministic competence state and the metastable state is smallest.

Equations

$$\begin{array}{rl}\frac{dk}{dt}& =a_k+\frac{b_k{k}^n}{k_0^n+k^n}-\frac{k}{1+k+s},\\ \frac{ds}{dt}& =\frac{b_s}{1+{(k/k_1)}^p}-\frac{s}{1+k+s}.\end{array}$$

$$\begin{array}{rl}dk& =f_1(k,s)dt+{ϵ}_{k}d{L}_t^1,\\ ds& =f_2(k,s)dt+{ϵ}_{s}d{L}_t^2,\end{array}$$

$${\nu }_{\alpha }(\mathrm{d}k)=C_{\alpha }|k{|}^{-(1+\alpha )}\mathrm{d}k,$$

$$\mathcal{A}u(k,s)$$

$$-u(k,s)][{ϵ}_k^{\alpha }{\nu }_{\alpha }(d{k}^{\prime }){\delta }_0(d{s}^{\prime })+{ϵ}_s^{\alpha }{\nu }_{\alpha }(d{s}^{\prime }){\delta }_0(d{k}^{\prime })]$$

$$=f_1{u}_k+{ϵ}_k^{\alpha }{\int }_{\mathbb{R}\setminus \{0\}}[u(k+k^{\prime },s)-u(k,s)]{\nu }_{\alpha }(d{k}^{\prime })$$

$$+f_2{u}_s+{ϵ}_s^{\alpha }{\int }_{\mathbb{R}\setminus \{0\}}[u(k,s+s^{\prime })-u(k,s)]{\nu }_{\alpha }(d{s}^{\prime }),$$

$$p_t(k,s,t)$$

$$-(f_{2}p{)}_s+{ϵ}_s^{\alpha }{\int }_{{\mathbb{R}}^1\setminus \{0\}}[p(k,s+s^{\prime },t)-p(k,s,t)]{\nu }_{\alpha }(d{s}^{\prime }),$$

$$p_t(k,s,t)$$

$$-\frac{{ϵ}_k^{\alpha }{C}_{\alpha }}{\alpha }\left[\frac{1}{(k-a{)}^{\alpha }}+\frac{1}{(b-k{)}^{\alpha }}\right]p(k,s,t)+{ϵ}_k^{\alpha }{C}_{\alpha }{\int }_{a-k}^{b-k}\frac{p(k+k^{\prime },s,t)-p(k,s,t)}{|k^{\prime }{|}^{1+\alpha }}\mathrm{d}{k}^{\prime }$$

$$-\frac{{ϵ}_s^{\alpha }{C}_{\alpha }}{\alpha }\left[\frac{1}{(s-c{)}^{\alpha }}+\frac{1}{(d-s{)}^{\alpha }}\right]p(k,s,t)+{ϵ}_s^{\alpha }{C}_{\alpha }{\int }_{c-y}^{d-y}\frac{p(k,s+s^{\prime },t)-p(k,s,t)}{|s^{\prime }{|}^{1+\alpha }}\mathrm{d}{s}^{\prime },$$

$$P_t(v,w,t)$$

$$-\frac{C_{\alpha }}{\alpha }(\frac{2{ϵ}_k}{b-a}{)}^{\alpha }\left[\frac{1}{(1+v{)}^{\alpha }}+\frac{1}{(1-v{)}^{\alpha }}\right]P+C_{\alpha }(\frac{2{ϵ}_k}{b-a}{)}^{\alpha }{\int }_{-1-v}^{1-v}\frac{P(v+v^{\prime },w,t)-P(v,w,t)}{|v^{\prime }{|}^{1+\alpha }}\mathrm{d}{v}^{\prime }$$

$$-\frac{C_{\alpha }}{\alpha }(\frac{2{ϵ}_s}{d-c}{)}^{\alpha }\left[\frac{1}{(1+w{)}^{\alpha }}+\frac{1}{(1-w{)}^{\alpha }}\right]P+C_{\alpha }(\frac{2{ϵ}_s}{d-c}{)}^{\alpha }{\int }_{-1-w}^{1-w}\frac{P(v,w+w^{\prime },t)-P(v,w,t)}{|w^{\prime }{|}^{1+\alpha }}\mathrm{d}{w}^{\prime }.$$

$$\begin{array}{rl}\frac{\mathrm{d}{P}_{i,j}}{\mathrm{d}t}& =-\frac{2}{b-a}[(f_{1}P{)}_{v,i}^{+}+(f_{1}P{)}_{v,i}^{-}]-\frac{2}{d-c}[(f_{2}P{)}_{w,j}^{+}+(f_{2}P{)}_{w,j}^{-}]\\ & -\frac{C_{\alpha }}{\alpha }(\frac{2{ϵ}_k}{b-a}{)}^{\alpha }\left[\frac{1}{(1+v_i{)}^{\alpha }}+\frac{1}{(1-v_i{)}^{\alpha }}\right]P_{i,j}-\frac{C_{\alpha }}{\alpha }(\frac{2{ϵ}_s}{d-c}{)}^{\alpha }\left[\frac{1}{(1+w_j{)}^{\alpha }}+\frac{1}{(1-w_j{)}^{\alpha }}\right]P_{i,j}\\ & +C_{hx}\frac{P_{i-1,j}-2P_{i,j}+P_{i+1,j}}{h^2}+C_{hy}\frac{P_{i,j-1}-2P_{i,j}+P_{i,j+1}}{h^2}\\ & +C_{\alpha }(\frac{2{ϵ}_k}{b-a}{)}^{\alpha }h\sum _{k_1=-I-i,k_1\ne 0}^{I-i}\frac{P_{i+k_1,j}-P_{i,j}}{|{v_k}_1{|}^{1+\alpha }}+C_{\alpha }(\frac{2{ϵ}_s}{d-c}{)}^{\alpha }h\sum _{k_2=-I-j,k_2\ne 0}^{I-j}\frac{P_{i,j+k_2}-P_{i,j}}{|{w_k}_2{|}^{1+\alpha }},\end{array}$$

$$\begin{array}{rl}{P}^{(1)}& =P^n+\Delta tR(P^n),\\ P^{(2)}& =\frac{3}{4}{P}^n+\frac{1}{4}{P}^{(1)}+\frac{1}{4}\Delta tR(P^{(1)}),\\ P^{n+1}& =\frac{1}{3}{P}^n+\frac{2}{3}{P}^{(2)}+\frac{2}{3}\Delta tR(P^{(2)}),\end{array}$$

Full text

label1label1footnotetext: This work is supported by NSFC (Grant Nos. 11531006, 11371367, 11271290, 11771062) and the Fundamental Research Funds for the Central Universities’, HUST: 0118011076.

1 Introduction

1.1 Background

A number of studies in different fields were concerned about the role of noise playing in dynamical systems, since noise may induce various delicate effects [1, 2, 3, 4, 5, 6]. Noise-induced transition phenomena were regarded in various dynamical systems including gene regulatory networks [7, 8], oscillators and electronic transport in physics [9, 10, 11, 12, 13], chemical reacting systems [14, 15], cancer cell proliferation and epidemiological models in pathology [16, 17].

These applications have been considered for Guassian noise. But when fluctuations are present in certain events, such as burst-like events, Gaussian noise is not appropriate. In this case, it is more appropriate to model the random fluctuations by a non-Gaussian Lévy motion with heavy tails and bursting sample paths. Non-Gaussion Lévy noise has found a lot of applications in different areas. Perc applied the non-Gaussian Lévy noise in economy and studied the impact of an important class of external perturbations on the evolution of cooperation in a spatially extended prisoner’s dilemma game [18, 19]. Lévy noise was also associated with saving energy in animals, e.g. albatross flight patterns, as reviewed in [20]. Zheng et. al considered how Lévy noise may affect gene regulatory networks in [21]. For further applications about non-Gaussion Lévy noise cf. [22, 23, 24, 25].

Gene expression is a noisy process [12, 26]. There are various studies about Gaussian noise in gene regulation [27, 28, 29, 30, 31]. However, the transcription of a gene can be a discontinuous or burst-like event, and mRNA is synthesized in intermittent but intense pulses or bursts [21, 22, 32, 33, 34]. Under this circumstance, Gaussian noise is not proper to model this phenomenon, but a stochastic process with discontinuous and heavy tail distribution, i.e., non-Gaussian Lévy noise [35, 36, 37, 38, 39, 40] appears more appropriate in describing the fluctuations in gene regulation.

In this study, we devote to studying the transitions in a two-dimensional genetic regulatory model driven by non-Gaussian $\alpha$-stable Lévy noises, which describe fluctuations with features such as heavy tails and jumps. The corresponding genetic regulatory model was established by Süel et al. [7]. Comparing with the experimental data, they demonstrated the reasonability and validity of the genetic regulatory model and considered a set of excursion trajectories in an excitable case under Gaussian noise. Instead of using the mean first exit time and first escape probability to study the bistable case of this model driven by $\alpha$-stable Lévy noise [41], we examine the most probable trajectories rather than stochastic trajectory sample paths to characterize the dynamical behaviors of the model. The nonlocal Fokker-Planck equation of stochastic systems with Lévy noise can be numerically simulated [42, 43, 44, 45, 46, 47, 48]. The noise induced most probable trajectories, which describe the MeKS network from the low ComK protein concentration (vegetative state) to the high ComK protein concentration (competence state), are obtained by examining the solution of the nonlocal Fokker-Planck equation. We use here numerical methods to solve the nonlocal Fokker-Plank equation, (for more details see Appendixes). We present choices of the non-Gaussianity index $\alpha$ and the noise intensity $\epsilon$ which induce the ComK protein transfer from the low concentration to the high concentration. Then we will analyse how different initial concentrations around the low ComK protein concentration affect the metastable state. Furthermore, we determine the optimal $\alpha$ and $\epsilon$ such that the distance between the deterministic competence state and the metastable state is smallest.

1.2 A stochastic model for the MeKS network

To investigate the dynamics of competence induction, Süel et al. [7] constructed a MeKS network:

[TABLE]

The symbols $k$ and $s$ represent the concentration of the ComK and ComS proteins, respectively. The parameters $a_{k}$ and $b_{k}$ stand for the basal and fully activated rates of ComK protein production, respectively. The parameter $k_{0}$ denote the concentration of the ComK protein needed for a $50\%$ activation. The Hill coefficients $n$ and $p$ are the cooperativities of the ComK protein auto-activation and ComS protein repression, respectively. The maximal expression rate of the ComS protein is $b_{s}$, and when $k=k_{1}$ the ComS protein attains its half-maximal rate. We choose suitable parameters [7]: $a_{k}\!=\!0.004,~{}b_{k}\!=\!0.14,~{}b_{s}\!=\!0.68,~{}k_{0}\!=\!0.2,~{}k_{l}\!=\!0.222,~{}n\!=\!2,~{}p\!=\!5$. In this case, the MeKS network has three equilibria: two stable equilibria (0.015262, 2.1574) (nodal sink corresponding to the low vegetative state), (0.15732,1.5781) (spiral sink corresponding to the high competence state) and one unstable equilibrium (0.08568,2.2469)(saddle point).

The core competence genetic circuit of the MeKS network is shown in FIG. 1. In B. subtilis population, mostly, the ComK protein is expressed by the comK gene at a relative low level. Thus, the ComK protein concentration is low, which corresponds to the vegetative state. Under circumstance of nutrient limitation, a small part of B. subtilis differentiate into the competence state, in which the ComK protein concentration is high. In this case, the cell can take in extracellular DNA from the environment [7, 49, 50, 51, 52]. It was found that the key transcription factor, the ComK protein, can activate the transcription of several genes which is necessary for the competence development. When the ComK protein concentration is high, then the competence develops [53, 54], from which we deduce that the ComK protein occupies a core position in the competence-signal-transduction network.

As discussed in the preceding subsection, gene expression in cells is subject to random fluctuations [12, 26], and such fluctuations are under real conditions typically non-Gaussian noise. Hence we consider the MeKS network under the influence of $\alpha$-stable Lévy motion, which can be written as the following stochastic differential equations:

[TABLE]

where $f_{1}(k,s)={{a}_{k}}+\tfrac{{{b}_{k}}{{k}^{n}}}{k_{0}^{n}+{{k}^{n}}}-\tfrac{k}{1+k+s},f_{2}(k,s)=\tfrac{{{b}_{s}}}{1+{{(k/{{k}_{1}})}^{p}}}-\tfrac{s}{1+k+s}.$ Here $(\epsilon_{k},\epsilon_{s})$ denote the diffusion coefficients, and $(L_{t}^{1},L_{t}^{2})$ is the two-dimensional Lévy motion with independent scalar symmetric $\alpha$-stable Lévy motions which have the same generating triplet $(0,0,\nu_{\alpha})$. Here the jump measure is a Borel measure in $\mathbb{R}\!\setminus\!\{0\}$ and defined by

[TABLE]

where $\displaystyle{C_{\alpha}=\frac{\alpha}{2^{1-\alpha}\sqrt{\pi}}\frac{\Gamma(\frac{1+\alpha}{2})}{\Gamma(1-\frac{\alpha}{2})}}$, $\alpha\in(0,2)$ is the non-Gaussianity index. For more information about $\alpha$-stable Lévy motions, see [2, 3].

In the following discussions, we use the scale transformation $k^{\prime}=10k,~{}s^{\prime}=2s$. For the noise intensity, we set $\epsilon=\epsilon_{k}=\epsilon_{s}$. After this transformation, the two stable equilibria become $(k_{-},s_{-})=(0.15262,4.3148)$ (the low concentration state) and $(k_{+},s_{+})=(1.5732,3.1562)$ (the high concentration state), and the unstable equilibrium becomes $(k_{u},s_{u})=(0.8568,4.4938)$. We define $D_{1}\!=\!(0,k_{u})$ as the low concentration region. If the most probable trajectories reach the high concentration region, where the ComK protein concentration is high, the competence develops.

This paper is organized as follows. In Section 2, we present the method how to compute the most probable trajectories. In Section 3, we discuss the most probable dynamics of the MeKS network with geometric tools (most probable trajectories) under stable Lévy noise. We finish this paper with conclusions and discussions in Section 4.

2 Method

When it comes to trajectories for a stochastic dynamical system, there is an apparent option of consideration: the almost sure trajectories [3, 55]. However, as we know, the almost sure trajectories, i.e. the sample trajectories (FIG. 2), which look like “noodles” in the phase plane, could hardly provide helpful information for understanding the system’s dynamics. Furthermore, to acquire comprehension of a stochastic dynamical phenomenon, stationary probability density functions for solution trajectories have been widely utilised [56, 57].

Notice that each sample trajectory is an “outcome” of a trajectory of a stochastic dynamical system, from a given initial state. Now, the question arises, which trajectory does the system most probably to evolve in phase plane? Our geometric tool, the most probable trajectories [3, 55], will tackle this problem. As is known, the solution of a stochastic system (1.2) is a stochastic process, and the probability density function of the stochastic process is governed by the corresponding Fokker-Planck equation [42, 57]. On the basis of it, we can get the most probable trajectory by computing the maximizer of the probability density function $p(k,s,t)$ at every time $t$. This offers geometric representations of most probable trajectories of a stochastic dynamic system [55].

First we derive the generator and the Fokker-Planck equation of system (1.2) (for more details see Section Symmetric $\alpha$-stable Lévy motions). Then we use the finite difference method to solve the Fokker-Planck equation (more details are given in Section Fokker-Planck equation). Plots of the so calculated probability density functions are shown in FIG. 3. We set the low concentration state $(k_{-},s_{-})$ as the initial point of probability density function. FIG. 3 shows that the peaks transfer from one stable point to another with increasing time. When time is small $(t=1)$, there is one peak around the low concentration. As time goes on, there are two peaks for the probability density, one peak is transferred to the competence region, where the ComK protein concentration is high, the other peak still stays around the vegetative region. For $t=20$, the probability density function has attained its stationary solution, finally it stays at the competence state (the high ComK concentration). For a fixed time $t_{i}$, we compute the maximum of $p(k,s,t_{i})$ to get the position $(k_{i},s_{i})$. We connect this series of $\{(k_{i},s_{i}),~{}i=1,~{}2,\cdot\cdot\cdot\}$, to get the most probable trajectory (from a given initial state). The time instants $t_{i}$’s need to be taken close enough, in order to get a reasonable approximation of the most probable trajectory.

Through computing the most probable trajectories, we analyze how the ComK protein transfers from the low concentration regime to the high one, with different noise parameters. Then we want to know when the most probable trajectories go over the saddle point. We define this time as the tipping time. Moreover, we will determine the choices of $\alpha$ and $\epsilon$ to benefit the development to competence.

3 Results

In this section, we analyse how the noise intensity $\epsilon$ and non-Gaussianity index $\alpha$ affect the most probable trajectories. As seen from FIG. 3, the probability density function attains its stationary state as time goes on. In the following, we choose $T=100$ as the computational terminal time.

Most probable trajectories and corresponding time series of the ComK protein for various noise intensity $\epsilon$ and non-Gaussianity index $\alpha$. The left plots of FIG. 4 show the impacts of $\epsilon$ and $\alpha$ on the most probable trajectories. Note that the red star is the low concentration state, the green star denotes the saddle point and the blue star is the high concentration state. The right plots of FIG. 4 show the impacts of $\epsilon$ and $\alpha$ on the corresponding time series of ComK and the red star is the high concentration state of ComK protein.

When $\epsilon$ is small, the most probable trajectories stay around the low concentration state (FIGs. 4 (a1) and (a2) ). Hence, a small noise intensity is not in favor of the expression of comK gene. When $\epsilon=0.15$ (see FIGs. 4 (b1) and (b2) ), the most probable trajectories first get near the unstable point, and then jump to the high concentration region. As time goes on, the most probable trajectories transit back to the low concentration region which means that the high concentration region is less stable or robust than the low concentration region in this case ($\epsilon=0.15$). With the increase of $\epsilon$, trajectories with $\alpha>0.25$ escape to the high ComK protein concentration region. But trajectories with $\alpha=0.25$ can not escape to the high concentration region by $T=100$. Note that when $\alpha=0.25$, the Lévy noise has larger jumps but lower jump frequencies. FIGs. 4 (c1) and (d1) exhibit that when $\epsilon=0.25$ and $0.4$, trajectories with $0.25<\alpha<2$ lead the ComK protein evolve to the high concentration region, and in this case, the jump frequencies are higher than for Lévy noise with $\alpha=0.25$.

Furthermore, we see that when the most probable trajectories get close to the saddle point

$(0.8568,4.4938)$, they will reach the high concentration soon. This is in accordance with the observation that the saddle point repulses the most probable trajectories, while the high concentration stable state appears to attract them. If the most probable trajectories pass the saddle point, then they arrive at the high concentration. We define the time spending from the low concentration state to the saddle point as the tipping time.

When the noise intensity is large, the tipping time becomes short (see FIGs. 4 (a2), (b2), (c2) and (d2) ). In the following, we will discuss the tipping time for different $\epsilon$ and $\alpha$.

Tipping time with different noise intensity $\epsilon$ and non-Gaussianity index $\alpha$. FIG. 5 shows the tipping time from the low concentration region to high concentration region with different $\epsilon$ and $\alpha$. It is worth pointing out that the most probable trajectories can not reach the high concentration region when the tipping time is up to 30.

As shown in FIG. 5 (a), when $\epsilon<0.2$, the tipping time first decreases and then increases as $\alpha$ increases. This means that $\alpha\approx 0.8$ with larger jumps is beneficial to making the expression of the comK gene. When $\epsilon>0.2$, the tipping time decreases as $\alpha$ increases, i.e. large $\alpha$ (smaller jumps with higher frequency) helps the expression of the comK gene. In this case, we may infer that the transition will occur at earlier times for larger $\alpha$. These results coincide with results in [35].

FIG. 5 (b) presents the tipping time as a function of $\epsilon$. The transition does not occur for $\alpha=0.25,$ due to its low jump frequency. For $\alpha>0.25,$ the tipping time decreases as the noise intensity increases. Except for $\alpha=0.25$, the tipping time decreases as $\alpha$ increases at a given noise intensity, because the jump frequency increases. In summary, a larger noise intensity with a larger non-Gaussianity index will play a more positive role in the expression of the comK gene and a relative small noise intensity with non-Gaussianity index closing to $1$ will play a positive role in the expression of the comK gene.

Combination of the noise intensity $\epsilon$ and the non-Gaussianity index $\alpha$ control the competence development. Now we discuss which combination of $\alpha$ and $\epsilon$ can make the ComK protein efficiently transit to the high concentration, i.e. the competence develops. It is shown in FIG. 6, that the region $L-L$ below the red curve can not make a ComK protein transit from the low concentration to the high one, but that the region $L-H$ above the red curve can induce the ComK protein transit to the high concentration. This figure provides clearly how to choose $\epsilon$ and $\alpha$ to develop the competence. Larger $\alpha$ with high frequency jumps is propitious to the expression of the comK gene. Small $\epsilon$ is not beneficial to the expression of the comK gene.

Various initial concentrations evolve to a metastable state. Here we study the effect of different initial concentrations near the nodal sink (the low ComK vegetative state) on the most probable trajectories (FIG. 7). As in FIG. 4, the red star denotes the nodal sink, the green star represents the saddle and the blue star is the spiral sink. As shown in FIGs. 7 (a1), (a2), when $\epsilon$ and $\alpha$ are fixed, we choose nine different initial concentrations surrounding the nodal sink, and find that the most probable trajectories show a slight tendency to the nodal sink, then ultimately gather at a metastable state, i.e. reaching the equilibrium in the sense of stochasticity. FIGs. 7 (b1), (b2) display the most probable time series of the ComK protein, from which we can see that the ComK protein tends to the metastable state as time goes on. Interestingly, as shown in FIG. 7 (b1), when $\epsilon$ is relatively small ($\epsilon=0.15$), the most probable trajectories of the nine different initial concentrations first attain the competence region (the high ComK concentration), then they return to the vegetative region (the low ComK concentration), and finally gather at a low metastable state. However, when $\epsilon$ becomes larger ($\epsilon=0.3$), the most probable trajectories of the nine different initial concentrations reach the competence region (the high ComK concentration), and gather at a relative high metastable state finally. The differences among the nine most probable trajectories become narrow initially, then spread relatively strong, and eventually they converge at the high metastable state.

The distance between the deterministic competence state and metastable state with different $\alpha$ and $\epsilon$. As shown in FIGs. 7 (a1) and (a2), the most probable trajectories of different initial concentrations gather at different metastable states under the influence of Lévy noise with different $\alpha$ and $\epsilon$. That is, the distance between the deterministic competence state and the metastable state varies with different $\alpha$ and $\epsilon$, from which we know that different Lévy noise induce various effects on the behavior of dynamical systems. Now we discuss which $\alpha$ and $\epsilon$ can make the distance smallest, i.e. acting most efficiently. We consider the distance between the deterministic competence state $(k_{+},s_{+})$ and the metastable state $(k(T),s(T))$ as $\textbf{d}=\sqrt{(k_{+}-k(T))^{2}+(s_{+}-s(T))^{2}}$. It is shown in FIG. 8, that d decreases firstly and then increases with the increase of $\alpha$ when the noise intensity is fixed. When the noise intensity becomes stronger, the minimum of d is bigger. When $\alpha=1.85$ and $\epsilon=0.2$, d is the smallest. Hence, Lévy noise with $\alpha=1.85$ and $\epsilon=0.2$ have the smallest effect on the deterministic competence state.

4 Conclusion

We have studied the most probable trajectories of a two-dimensional system (MeKS network) under the influence of non-Gaussian stable Lévy noise. On the basis of the nonlocal Fokker-Planck equation, we have described the most probable trajectories of the MeKS network from the low ComK protein concentration (vegetative state) to the high ComK protein concentration (competence state) under stable Lévy noise. It has been recently demonstrated that the non-Gaussian stable Lévy noise is more appropriate to describe the burst-like phenomenon in synthesis process of biochemical systems. Therefore, we here study the most probable trajectories of the MeKS network under the influence of stable Lévy noise.

The most probable trajectories of the MeKS network influenced by Lévy noise with various noise intensity $\epsilon$ and non-Gaussianity index $\alpha$ have been investigated. We have found that for fixed $\alpha>0.25,$ the ComK protein stays at the low concentration region for small noise intensity, and it begins to transit to the high concentration region with the increase of $\epsilon$. Also, the tipping time from the low concentration region to the high concentration region as a function of $\alpha$ with various $\epsilon$ has been analysed. We have discovered that for fixed $\epsilon$, the tipping time firstly decreases and then increases as $\alpha$ increases. Thus, we can determine an optimal combination of $\alpha$ and $\epsilon$ to make the tipping time shortest. In addition, combinations of parameters of $\alpha$ and $\epsilon$ which generate a ComK protein transit from the low concentration to the high concentration have also been revealed. Furthermore, we have selected different initial concentrations around the low ComK protein concentration to examine the most probable trajectories of the MeKS network, to discover that the most probable trajectories gather at a metastable state. Lévy motion with a small $\alpha$ ($\alpha\leq 0.25$) which has a lower jump frequency can not make the ComK protein transit to the competence region. For small $\epsilon$, Lévy motion with a relative small $\alpha$, which has larger jumps, is more convenient to make the ComK protein transit to the competence region. However, for larger $\epsilon$, Lévy motion with a larger $\alpha$ which has more small jumps is more convenient to transit to the competence region. Additionally, when $\alpha=1.85$ and $\epsilon=0.2$, the distance between the deterministic competence state and metastable state is smallest.

In short, we have utilised this geometric tool, the most probable trajectories, to visualise the dynamics of the MeKS networks, which offer an optimal path for the ComK protein escapes from the low concentration (vegetative region) to the high concentration (competence region).

Appendixes

Symmetric $\alpha$-stable Lévy motions

A stochastic process $\{L_{t}:t\geq 0\}$ defined on $\mathbb{R}^{d}$ is a symmetric $\alpha$-stable Lévy process with $0<\alpha<2$ if the following conditions are satisfied:

1. (1)

$L_{0}=0$, a.s.; 2. (2)

For any choice of $n\geq 1$ and $0<t_{0}<t_{1}<\cdot\cdot\cdot<t_{n}$, random variables $L_{t_{0}},~{}L_{t_{1}}-L_{t_{0}},~{}L_{t_{2}}-L_{t_{1}},\cdot\cdot\cdot,~{}L_{t_{n}}-L_{t_{n\!-\!1}}$ are independent (independent increments property); 3. (3)

The distribution of $L_{s+t}-L_{s}$ does not depend on $s$ (temporal homogeneity or stationary increments property); 4. (4)

$L_{t}$ has stochastically continuous sample paths, i.e., for every $s>0,~{}L_{t}\rightarrow L_{s}$ in probability, as $t\rightarrow s$.

Fokker-Planck equation

The generator $\mathscr{A}$ of the process $(k_{t},s_{t})$ is defined as (see [2, 23, 41]),

[TABLE]

where $\delta_{0}$ is the delta measure concentrated at [math] with property $\int_{-\infty}^{\infty}g(z)\delta_{0}(dz)=g(0)$.

The Fokker-Planck equation for the distribution of the conditional probability density

$p(k,s,t)=p(k,s,t|k_{0},s_{0},0)$, i.e., the probability of the process $(k_{t},s_{t})$ has value $(k,s)$ at time $t$ given it had value $(k_{0},s_{0})$ at time [math], is $p_{t}=\mathscr{A}^{*}p$ given by [1, 3], where $\mathscr{A}^{*}$ is the adjoint operator for $\mathscr{A}$.

The adjoint operator for $\mathscr{A}$ can be found out, as the integral part is a symmetric operator. Thus we have the Fokker-Planck equation:

[TABLE]

with the initial condition: $p(k,s,0)=\delta(k-k_{0},s-s_{0})$.

Numerical method

In the following, we present the numerical algorithms the case of $D=(a,b)\times(c,d)$ with the absorbing condition. Noting that the absorbing condition dictates that $p$ vanishes outside $D$, we can simplify Eq. (Fokker-Planck equation) as

[TABLE]

for $(k,s)\in D$; and $p(k,s,t)=0$ for $(k,s)\notin D$.

Set $v=\frac{2}{b\!-\!a}(k\!-\!a)\!-\!1$, $w=\frac{2}{d\!-\!c}(s\!-\!c)\!-\!1$, $p(\frac{b\!-\!a}{2}v\!-\!\frac{a\!-\!b}{2},\frac{d\!-\!c}{2}w\!-\!\frac{c\!-\!d}{2},t)=P(v,w,t)$, i.e., $k=\frac{b\!-\!a}{2}v\!-\!\frac{a\!-\!b}{2}$, $s=\frac{d\!-\!c}{2}w\!-\!\frac{c\!-\!d}{2}$, then we get

[TABLE]

for $(v,w)\in D^{\prime}=(-1,1)\times(-1,1)$; and $P(v,w,t)=0$ for $(v,w)\notin D^{\prime}$.

Then, we use a numerical method to discretize the nonlocal Fokker-Planck equation (4.4). For the spatial direction, the advection term $-(f_{1}P)_{x}$ and $-(f_{2}P)_{y}$ are discretized by the third-order WENO method given in [58] and the singular integral term are used a modified trapezoidal rule to approximate[42, 43]. We divide the interval $[-2,2]\times[-2,2]$ in space into $(4I)^{2}$ sub-intervals and define $v_{i}=ih$, $w_{j}=jh$ for $-2I\leq i,j\leq 2I$, where $h=1/I$. Denoting the numerical solution of $P$ at $(v_{i},w_{j},t)$ by $P_{i,j}$, we obtain the semi-discrete equation

[TABLE]

for $0<\alpha<2$ and $i,j=-I+1,\cdots,-2,-1,0,1,2,\cdots,I-1$, where the constant $C_{hx}\!=\!-C_{\alpha}(\frac{2\epsilon_{k}}{b\!-\!a})^{\alpha}\zeta(\alpha-1)h^{2-\alpha}$ and $C_{hy}\!=\!-C_{\alpha}(\frac{2\epsilon_{s}}{d\!-\!c})^{\alpha}\zeta(\alpha-1)h^{2-\alpha}$, $\zeta$ is the Riemann zeta function, and the $\pm$ superscripts denote the global Lax-Friedrichs flux splitting defined as $(f_{k}P)^{\pm}_{i,j}=\frac{1}{2}(f_{k}P_{i,j}\pm\alpha_{k}P_{i,j})$ with $\alpha_{k}=\max|f_{k}(x,y)|$, $k=1,2$. The absorbing condition requires that the values of $P_{i,j}=0$ if the index $|i|\geq I$ or $|j|\geq I$.

For time discretization, we use a third-order total variation diminishing Runge-Kutta method provided in [59]. In particular, for the ordinary differential equation $\dfrac{{\rm d}P}{{\rm d}t}=R(P)$, the method can be written as

[TABLE]

where $P^{n}$ denotes the numerical solution of $P$ at time $t=t_{n}$.

Acknowledgements

We would like to thank Dongfang Li, Ke Yin, Yayun Zheng, Xiujun Cheng and Rui Cai for helpful discussions.