L\'{e}vy noise-induced transitions in gene regulatory networks

TL;DR

This paper investigates how Lévy noise influences gene expression transitions in a two-dimensional B. subtilis model, revealing optimal noise conditions and the impact on basin stability of the competence state.

Contribution

It introduces a two-dimensional gene network model incorporating Lévy and Brownian noise, analyzing transition dynamics with new quantitative measures and control strategies.

Findings

Optimal noise conditions enhance transition probability.

Lévy noise significantly affects basin stability.

Transition metrics depend on noise strength ratio.

Abstract

Important effects of noise on a one-dimensional gene expression model involving a single gene have recently been discussed. However, few works have been devoted to the transition in two-dimensional models which include the interaction of genes. Therefore, we investigate here, a quantitative two-dimensional model (MeKS network) of gene expression dynamics describing the competence development in the B. subtilis under the influence of L\'evy as well as Brownian motions, where noises can do the B. subtilis a favor in nutrient depletion. To analyze the transitions between the vegetative and the competence regions therein, two deterministic quantities, the mean first exit time (MFET) and the first escape probability (FEP) from a microscopic perspective, as well as their averaged versions from a macroscopic perspective, are applied. The relative contribution factor (RCF), the ratio of non-Gaussian and Gaussian noise strengths, is adopted to implement optimal control in these transitions. Schematic representations indicate that there exists an optimum choice that makes the transition occurring at the highest probability. Additionally, we use a geometric concept, the stochastic basin of attraction, to exhibit a pictorial comprehension about the influence of the L\'{e}vy motion on the basin stability of the competence state.

Equations

$$\displaystyle\hskip 65.44142pt\left\{\begin{array}[]{ll}\Psi=\varepsilon_{k}=\varepsilon_{s},\sigma_{k}=\sigma_{s}=0.&\hbox{$\alpha\in(0,2)\leavevmode\nobreak\ \leavevmode\nobreak\ (\alpha-stable\leavevmode\nobreak\ L\acute{e}vy\leavevmode\nobreak\ motion)$,}\\ \Psi=\sigma_{k}=\sigma_{s},\varepsilon_{k}=\varepsilon_{s}=0.&\hbox{$\alpha=2\leavevmode\nobreak\ \leavevmode\nobreak\ (Brownian\leavevmode\nobreak\ motion)$.}\end{array}\right.$$

$$\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}\frac{dk}{dt}& =f(k,s)+{\sigma }_k{\dot{B}}_t^1+{\epsilon }_k{\dot{L}}_t^1,\\ \frac{ds}{dt}& =g(k,s)+{\sigma }_s{\dot{B}}_t^2+{\epsilon }_s{\dot{L}}_t^2.\end{array}$$

$$\begin{array}{rl}\mathcal{A}u(k,s)& =-1,(k,s)\in D,\\ u(k,s)& =0\text{leavevmode}\text{leavevmode}\text{leavevmode}(k,s)\in D^c,\end{array}$$

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

$$\begin{split}\hskip 142.26378pt\mathscr{A}p(k,s)&=0,\qquad(k,s)\in D,\\ p(k,s)&=\left\{\begin{array}[]{ll}1,&\hbox{$(k,s)\in E$;}\\ 0,&\hbox{$(k,s)\in D^{c}\setminus E$.}\end{array}\right.\end{split}$$

$$\mathcal{A}v(z,w)$$

Full text

Introduction

Random fluctuations are ubiquitous in real life, and may change the fundamental dynamical properties of a system. A number of studies have reviewed that noises have rather counterintuitive and constructive effects on the behavior of a dynamical system. Numerous research in different fields has demonstrated that noise can play a pivotal role in dynamical systems in biology, physics, geophysics, chemistry, finance or engineering. This topic has attracted more and more attention in recent years [4, 5, 6, 1, 3, 2]. The present contribution is devoted to investigating the influence of a combination of Gaussian and non-Gaussian noises on a gene regulatory network, where the interaction of genes is included in substantial extensions to recent work.

It is known that messenger RNA (mRNA) and protein molecules are synthesized during the process of transcription and translation, respectively, which constitutes major procedures of gene expression. Bioscientists discovered that key transcription and regulatory factors pulse on and off repeatedly, often stochastically, even when cells are maintained in constant conditions. This type of spontaneous dynamic behavior is pervasive, appearing in diverse cell types from microbes to mammalian cells. The initiation of transcription, i.e., the binding of RNA polymerase at the promoter region of DNA, is inherently stochastic in biochemical reactions. Furthermore, the variability in the number of regulatory molecules provides noise sources from the extrinsic environment. Together, the inherent uncertainty in gene expression dynamics and external noisy environment make up the intrinsic and extrinsic noise sources [7, 8, 11, 9, 10]. Given that, it is necessary to take into account stochastic fluctuations when considering the dynamics of a gene expression system. It is often assumed that the noise is Gaussian. This arises due to the assumption that the perturbation is the result of a large number of independent interactions of bounded strength. However, this assumption is often not appropriate, in particular when the fluctuations are abrupt pulses or extreme events. Just as Levine et al [12] pointed out that, pulsatile dynamics can produce “long-tailed” distributions in static measurements based on flow cytometry and microscopy snapshots. Pulsatile dynamics is emerging as a central, and still largely unexplored area in the cell. In this situation, it is more appropriate to model the fluctuations by a process with heavy tails and discontinuous sample paths [5, 6, 13]. Lévy motion is an appropriate prototype for non-Gaussian processes with jumps. Therefore, we will here study such a basic gene regulatory network under the influence of Lévy motion.

Most of the B. subtilis population is in the vegetative state in which the comK gene expresses ComK protein at a low level. Under nutrient limitation, only a small fraction of B. subtilis transiently differentiates into the competence state in which the ComK protein at a high level, permitting the cell to take up extracellular DNA and to incorporate it into their own chromosome [14, 15, 16, 17], which is regulated by the MeKS network [18]. The core components of the MeKS network are shown in Fig.1a: the autoregulatory positive feedback loop regulated by the ComK protein, the ComS protein’s competitive interference with the proteolytic degradation of the ComK protein, and the indirect repression of the comS gene by the ComK protein [19, 20, 18]. The ComK protein can activate its own expression, shaping a positive feedback loop [21, 22]. Sinderen [23] proposed that, the key transcription factor ComK protein, functioning as a master regulator, can activate the transcription of several genes which is necessary for the competence. When the concentration of the ComK protein is high, the competence develops [25, 24], from which we know that the ComK protein occupies a central position in the competence-signal-transduction network. The relationship between the ComK and the ComS protein is responsible for the formation of a negative feedback loop. Together, positive and negative feedback loops interacting with each other constitute the MeKS network.

In order to examine the dynamics of competence development governed by the MeKS network structure, basing on experimental observations, Süel et al [18] built a mathematical model (2) given in Methods, which involves both the direct positive and the ComS-related negative feedback loops of ComK. In the MeKS model (2), the symbols $k$ and $s$ denote the concentration levels of the ComK and ComS proteins, respectively. The basal and fully activated rates of ComK production are $a_{k}$ and $b_{k}$, respectively. The parameter $k_{0}$ represents the concentration of the ComK protein needed for a $50\%$ activation. The cooperativities of ComK auto-activation and ComS repression are indicated by the Hill coefficients $n$ and $p$, respectively. The ComS protein has the maximum expression rate $b_{s}$ and the half-maximal rate at $k=k_{1}$. The ComS protein competitively inhibits the degradation activity of the MecA-ClpP-ClpC complex that targets the ComK protein, which is responsible for the nonlinear degradation terms described at the end of the MeKS model (2).

To understand the dynamics of the MeKS network without noise, with parameters suggested in ref.18, we calculate the nullclines in the ComK-ComS phase plane. Fig.1b exhibits three equilibria: two stable equilibria (nodal sink and spiral sink), whose ComK concentrations are low (vegetative state) and high (competence state), and an unstable saddle at a medium-concentration of ComK. The phase diagram of the MeKS model in the $b_{k}$-$b_{s}$ plane is depicted in refs.9 and 18 , which delineates four regions: (i) a monostability region in which only one stable equilibrium exists. (ii) a bistability region with three equilibria, two of them are stable. (iii) an excitability region with three equilibria ,only one of them is stable. (iv) a limit cycle oscillations region with only one unstable equilibrium. Our focus, however, is on the problem for the transition between both stable states. Therefore, the bistability region (ii) is chosen for our study.

It was found that stochasticity in gene expression can induce random transitions between both stable steady states, giving rise to binary distributions in the cell population, i.e., two pronounced subpopulations with different protein concentrations appear, respectively [7]. Indeed, the influence of noise can cause random perturbations to the stable steady state of a system. Fluctuations sufficiently strong enable the system to escape it, residing in another stable steady state [9, 10]. The existence of two major ingredients, the positive feedback loop involved in dynamical system and the additional nonlinearity, contribute to the phenomenon of bistability [26, 27, 25].

According to previous works [9, 12, 28, 29, 30], noise is pervasive in the gene expression events. On occasion, these events show burst-like dynamical behaviors. Thus, Brownian motion and $\alpha$-stable Lévy motion is suitable to model the noise sources in gene regulatory systems. Here, the $\alpha$-stable Lévy motion has jumps, which is a special kind of non-Gaussian process defined by stable Lévy random variables. It has a larger jump magnitude with lower jump frequency when $\alpha$ closes to [math], while smaller jump magnitude with higher jump frequency when $\alpha\in(1,2).$ It is substantially different from Brownian motion whose paths are continuous but almost surely nowhere differentiable. The $\alpha$-stable Lévy motion’s paths are continuous from the right and have left limits (càdlàg) at each time [5, 6]. Notice that a càdlàg function can only have a countable number of jumps at most, permitting the system to transit between one stable state and another one which is enclosed in a non-adjacent domain.

To examine the noise-induced transition between both stable states, two deterministic quantities, the mean first exit time (MFET) and the first escape probability (FEP) from a microscopic perspective, as well as the average first exit time (AMFET) and the average first escape probability (AFEP) from a macroscopic perspective, are utilised here. They are governed by differential-integral equations (4) and (6) (see Methods), and a visualized representation for the MFET and the FEP in the two-dimensional stochastic MeKS model (3) are exhibited in Fig.2. The microscopic perspective indicates that we analyze the dynamical behavior of a system from a certain point $(k,s)$, but in the macroscopic perspective the entirety dynamical behavior in a domain $D$ is regarded. Effects of Lévy and Brownian motions on a one-dimensional gene expression model have recently been investigated [29]. Here, a more general two-dimensional MeKS model describing the competence development in the B. subtilis under the influence of Lévy as well as Brownian motion is studied. However, it becomes more challenging but difficult to conduct such a simulation in the two-dimensional model with Lévy motions. Nevertheless, we manage to calculate these deterministic quantities by accurate and efficient numerical methods [31, 13, 32, 34, 33, 35], instead of Monte-Carlo simulation.

The effects of noise on the model of autoregulatory gene expression involving a single gene has been discussed in ref.28. However, the aforementioned work is only for a single gene activating a positive feedback loop, but avoids any consideration with other genes. However, the interplay of different genes is crucial for real genomic systems. Pal et al [9] studied a one-dimensional and a two-dimensional model in gene expression dynamics in B. subtilis, respectively. They obtained some early signatures of regime shifts in gene expression dynamics for the one-dimensional case, such as critical slowing down (CSD), rising variance and a ratio of two mean first passage times, while for the two-dimensional model, they did not get as many results as for the one-dimensional case. So far, little work has been devoted to effects of noise on the transition for the two-dimensional model of the MeKS network in B. subtilis. Our main goal is to utilise the knowledge of stochastic dynamics [6] to gain some deeper insight into effects of various noises on the MeKS network, and in particular we try to explore and offer clues on the connection between vegetative growth and competence in the MeKS network.

Furthermore, it has been recently discovered that any disease progression can be divided into a normal state, a pre-disease state, and a disease state [37, 36], similar to the low vegetative state, the medium concentration state, and the high competence state in our work. As we know, diverse diseases pose a tremendous threat to people’s health, and it has become a major challenge in medicine and biology for cancer research. In complex diseases networks, Huang [38] pointed out the gene expression noise induces effects on the choices between the normal and the malignant phenotype of the cell during cancer progression. Also, Jia et al [39] investigated the stochastic phenotype switching within the bacterial populations, and proposed the importance of the identification of the critical state for an early diagnosis in complex diseases networks. However, they just considered the influence of the Gaussian noise. Here, through our exploration for transitions between both stable states in gene expression, we expect our work to provide clues for the prevention and therapy of complex diseases system with a normal state, a pre-disease state, and a disease state.

Results

To illuminate our work, we choose parameters of the bistability region (ii) as listed in ref.18 for our numerical experiments, namely, $a_{k}=0.004,b_{k}=0.14,b_{s}=0.68,k_{0}=0.2,k_{1}=0.222,n=2,p=5.$ Then, the deterministic MeKS model (2) has two stable states at $(k_{l},s_{l})\approx(0.015262,2.1574)$ (the low ComK vegetative state), $(k_{h},s_{h})\approx(0.15732,1.5781)$ (the high ComK competence state) and the unstable state at $(k_{u},s_{u})\approx(0.08568,2.2469)$ (the medium ComK state). The effects of noise on the one-dimensional model involving a single gene has been investigated [29, 28], and it has been recently found that the noise intensity and the size of jumps related with the Lévy motion index $\alpha$ can take effect on the concentrations within a cell in gene expression [29]. However, little work has been committed to the transition for the two-dimensional model. Here, we study a more general situation, i.e., a quantitative two-dimensional model (MeKS network) of gene expression dynamics describing competence development in B. subtilis with noises. Next, we present our results on the effects of Gaussian Brownian motion vs. non-Gaussian $\alpha$-stable Lévy motion on the first exit time and the first escape probability (both from microscopic and macroscopic perspectives) as well as basin stability for the competence state in the two-dimensional MeKS model.

Shorter MFET for higher noise intensity and smaller jump magnitude. For consideration of the residence time of the vegetative stable state, we fix the low vegetative region $D\!=\!(0,0.793)\!\times\!(0,2.46)$ enclosing the low ComK stable state $(k_{l},s_{l})$ from the ComK-ComS phase plane. Now, from microscopic and macroscopic perspectives respectively, we will investigate the influence of Brownian motion and $\alpha$-stable Lévy motion on the first exit time, which is the mean time of the ComK protein residing in the low vegetative region $D$ before exiting to the high competence region. In Methods, the differential-integral equations (4) offers how to get the mean first exit time $u(k,s)$, which is subject to the nonlocal “Dirichlet” condition, i.e. $u(k,s)$ vanishes outside of $D$.

Without stochastic perturbations, trajectories starting from the initial state around one stable state can not switch to the region enclosing another stable state. Thus, the residence time staying in the initial region is infinite. However, due to the influence of noise, the system can perform a given transition between both stable states.

Here, we examine the impact of the Brownian motion and $\alpha$-stable Lévy motion on the mean first exit time, separately. For convenience, we introduce a notation $\Psi$ as the noise intensity under different stochastic fluctuations, i.e.,

[TABLE]

First, we analyse the noise intensity and Lévy motion parameter from a microscopic viewpoint (MFET for the low ComK vegetative state). It is shown in Fig.3a that, when the noise intensity is small enough, MFET increases as $\alpha$ increases, and grows rapidly at $\alpha=2$ (special case corresponding to Brownian motion). That is, if the noise intensity is slight, then the vegetative state is more stable under Brownian motion than $\alpha$-stable Lévy motion. Nonetheless, when the noise intensity becomes larger, the tendency of MFET changes from increase to decrease as $\alpha$ increases, implying that there are critical values of $\alpha$ resulting in a change of MFET from increase to decrease. Meanwhile, when $\alpha$ is fixed at a certain value, there is a decline in MFET as the noise intensity $\Psi$ increases. In other words, $\Psi$ weakens the stability of the vegetative state. The macroscopic situation (AMFET for the low ComK vegetative region) is similar to the consequence in the microscopic situation, we omit its figure here. As shown in Fig.3a, a larger $\Psi$ induces a shorter residence time. In addition, compared with Brownian motion, the $\alpha$-stable Lévy motion makes MFET shorter when $\Psi$ are the same. Together, both microscopic and macroscopic situations demonstrate the consistency of the stability behaviour of the state $(k_{l},s_{l}).$

The pictorial representation of the aforementioned results enlighten us on the transition subject, i.e., if we expect the ComK protein to exit from the low vegetative state to the high competence state, then a higher noise intensity and a larger $\alpha$ (smaller jump magnitude with higher frequency) should be selected.

Optimal control on the stability of the vegetative state. These results point out different effects of Gaussian noise (Brownian motion) and non-Gaussian noise ($\alpha$-stable Lévy motion) on the first exit time, individually. Now, we would like to discuss the influence of a combination of both types of noise on the first exit time. Due to the vital role of the ComK protein in the competence development, here, we apply the relative contributor factor (RCF) $\lambda$ proposed in ref.29. In detail, $\lambda=\frac{\varepsilon_{k}}{\sigma_{k}}$, where $\sigma_{k}$ belongs to $(0,1]$, $\varepsilon_{k}$ stands for the non-Gaussian noise intensity and $\sigma_{k}$ for the Gaussian noise intensity, with the added constraints $\varepsilon_{k}=\varepsilon_{s},\sigma_{k}=\sigma_{s}$ and $\varepsilon_{k}+\sigma_{k}=1$.

Fig. 3b,c exhibit MFET around the low vegetative state (nodal sink) as a function of $\lambda$ for $\alpha=0.5$ and $\alpha=1.5$, respectively. From a microscopic viewpoint, we choose four different initial concentrations around the nodal sink, i.e., the red point is in close proximity to the nodal sink, the blue lies in left direction of the nodal sink, and the green and cyan in right direction of the nodal sink. As shown in Fig.3b, with the increase of $\lambda$, all initial points of MFET increase fast, after reaching their maxima, the MFET tends to a constant. Meanwhile, the critical point $\lambda_{0}$ corresponding to the maximum of MFET is much larger than $1$, implying that the Lévy motion makes a much stronger contribution on MFET than the Brownian one. However, from Fig.3c, when $\alpha=1.5,$ the critical point $\lambda_{0}$ shifts towards a much smaller value than for $\alpha=0.5,$ signifying that the Brownian motion takes a more significant effect on the MFET gradually. Note that there is an obvious decline in MFET before arrival at a constant concentration. Furthermore, compared with Fig.3b,c, MFET for $\alpha=0.5$ is much higher for $\alpha=1.5.$

Next, we turn to the macroscopic viewpoint. The graphic representation of Fig.3d is more or less the same as Fig.3b. It is worth noticing that, when $\lambda$ is fixed, AMFET takes on a rising trend as $\alpha$ decreases, which implies that a smaller $\alpha$ can stabilize the vegetative state. Simultaneously, due to the gradually significant effect on MFET by Brownian motion, the critical point $\lambda_{0}$ shifts towards a much smaller value as $\alpha$ increases.

Together, the critical point $\lambda_{0}$ corresponding to the maximum of MFET provides us the optimal control on the stability for the vegetative state. When the ratio of both noise intensities equals $\lambda_{0}$, it is more difficult to exit the vegetative region. In other words, noises determined by this ratio can enhance stability of the vegetative state. Actually, this is the phenomenon of noise- enhanced stability, a resonance-like effect, which has been found experimentally and theoretically in various systems[40, 41, 42, 43]. However, if we want to weaken the stability of the vegetative state, and expect the concentration of the ComK protein to acquire a high level, then a larger value of $\alpha$ is useful.

First escape probability.

Now, we employ another concept to discuss the impacts of the noise intensity and $\alpha$ on the behavior of the transition probability, i.e., the likelihood of the first switching to a escape-target region. Here, we are concerned with the escape probability between the low vegetative region and the high competence region. According to the escaping phenomenon, there are two scenarios of the escape-target region: one is adjacent to the boundary of the starting region, but the other is non-adjacent. The first escape probability governed by the differential-integral equations (6) can be solved by a numerical scheme efficiently in Methods.

Escape probability from the vegetative region to the adjacent competence region.

In this part, we are concerned with the escape probability from the vegetative region $D\!=\!(0,0.0793)\!\times\!(0,2.46)$ to the adjacent competence region $E\!=\!(0.0793,\infty)\!\times\!(0,2.46)$. From a microscopic perspective, Fig.4a,b,c depict FEP as a function of $\alpha$ with different noise intensities. Fig.4a describes the behavior of the escape probability of the left point of the nodal sink. As shown in the figure, for different noise intensities, as $\alpha$ increases, there is an increasing tendency in FEP initially, but after reaching its maximum, FEP declines slowly. Notice that FEP is smallest at $\alpha=2$ (Brownian motion). From Fig.4b, we find that the behavior of the FEP of the point near the nodal sink is similar, except the much slower decline rate. However, for the right point of the nodal sink (Fig.4c), the FEP increases as $\alpha$ increases. Interestingly, due to its relative high concentration, the smallest noise intensity ($\Psi=0.0125$) makes FEP decreasing at $\alpha=2$, however, other noise intensities lead FEP to increase. Compared with Fig.4a,b,c, we observe that the transition probability varies with different initial concentrations, and it increases when the ComK protein initial concentration $(k,s)$ is close to the relative high concentration.

Now, we turn to look at the whole escape behavior. It is shown in Fig.4d that, for the smallest noise intensity ($\Psi=0.0125$), the escape behavior is about the same as the point near the nodal sink in Fig.4b. While for other noise intensities, AFEP increases as $\alpha$ increases, the AFEP drops at $\alpha=2$ (Brownian motion).

Escape probability from the vegetative region to the non-adjacent competence region.

The escape probability from the vegetative region $D\!=\!(0,0.04)\!\times\!(1.26,2.46)$ to the non-adjacent competence region $E\!=\!(0.138,0.178)\!\times\!(1.26,2.46)$ including the spiral sink is considered in this part, and the areas of the two regions are same. There exists a critical value of $\alpha_{0}$ with the highest escape probability in Fig.4e,f. First, from a microscopic perspective, Fig.4e depicts FEP as a function of $\alpha$ with different noise intensities. As $\alpha$ increases, there is an increasing tendency in FEP initially, but after reaching its maximum, FEP drops down to zero. Notice that FEP is much smaller than the adjacent scenario, and it is zero at $\alpha=2$, due to the fact that Brownian motion’s trajectories are continuous, thus they can not jump. This last result is in accordance with the existing work[44, 45, 29]. Macroscopically, as shown in Fig.4f, the escape behavior of the whole vegetative region enclosing the nodal sink is similar.

Optimal control for FEP from the vegetative region to the adjacent competence region.

To implement the optimal control for the first escape probability of the transition from the vegetative region $D\!=\!(0,0.0793)\!\times\!(0,2.46)$ to the adjacent competence region $E\!=\!(0.0793,\infty)\!\times\!(0,2.46)$, also, the relative contributor factor $\lambda$ is utilised. Microscopically, in Fig.5a,b, the same four initial concentrations around the low vegetative state(nodal sink) in Fig.3b,c are selected. As shown in Fig.5a, when $\alpha=0.5$, before entering to their constant concentrations, there is an opposite tendency for different points with the increase of $\lambda$, i.e., FEP is decreasing for the two right points of the nodal sink, while the point near the nodal sink and the left point of the nodal sink increase. Meanwhile, the critical value $\lambda_{0}$ corresponding to minima or maxima is large for $\alpha=0.5$, which indicates that the Lévy motion brings much more contribution to the transition. Then, when $\alpha=1.5$, there is a similar tendency in Fig.5b with Fig.5a, except for the position of the critical value $\lambda_{0}$. Now, macroscopically, in Fig.5c, with the increase of $\lambda$, a slight increase in AFEP emerges for a much smaller $\lambda$, and after reaching its maximum, AFEP declines to arrive at a constant. When $\lambda$ is fixed at some value, AFEP increases as $\alpha$ increases. This indicates that we can choose a larger $\alpha$ (a smaller jump magnitude with a higher frequency) to make the transition occur at the highest probability.

Optimal control for FEP from the vegetative region to the non-adjacent competence region.

Following up with the optimal control for the adjacent scenario, the relative contributor factor $\lambda$ is utilised too, for the transition from the vegetative region $D\!=\!(0,0.04)\!\times\!(1.26,2.46)$ to the non-adjacent competence region $E\!=\!(0.138,0.178)\!\times\!(1.26,2.46)$. Microscopically, four initial concentrations around the low vegetative state (nodal sink) are selected. It is shown in Fig.5d, when $\alpha=0.5$, that all the selected points increase with the increase of $\lambda$. Meanwhile, the critical value $\lambda_{0}$ corresponding to the maximum is large for $\alpha=0.5$, implying that the Lévy motion gives much more contribution to the transition. When $\alpha=1.5$, the tendency in Fig.5e is the same with Fig.5d, but the position of the critical value $\lambda_{0}$ shifts towards a much smaller point. Macroscopically, as shown in figure Fig.5f, AFEP increases as $\lambda$ increases. When $\lambda$ is fixed at a certain value, AFEP increases as $\alpha$ decreases. Hence, the results indicate that $\alpha_{0}\approx 0.5$ (a larger jump magnitude with a lower frequency) is an ideal choice to make the transition to a high level at the highest probability.

FEP between the vegetative region and the competence region for different scenarios.

The preceding results exhibit the behavior of escape probability in one direction, i.e., from the vegetative region to the competence region. Süel et al [18] pointed that individual cells in B. subtilis can switch to an alternative state and, after some time, switch back again. It is necessary to survey the behavior of escape probability in both directions. We take notations FEP+ as the FEP from the vegetative state (nodal sink) to the competence region, and FEP- the FEP from the competence state (spiral sink) to the vegetative region. Fig.6 describes the behavior of the FEP in both directions for different scenarios from a microscopic perspective. Two regions with equal areas, the low vegetative region $D_{1}\!=\!(0,0.086291)\!\times\!(1.26,2.46)$ and the adjacent high competence region $\tilde{D}_{1}\!=\!(0.086291,0.172582)\!\times(1.26,2.46)\!,$ are selected in Fig.6a,b. And for the non-adjacent scenario (Fig.6c,d), we choose regions $D_{2}\!=\!(0,0.04)\!\times\!(1.26,2.46)$ and $\tilde{D}_{2}\!=\!(0.138,0.178)\!\times\!(1.26,2.46)$.

Fig.6a shows the behaviors of FEP in both directions when $\alpha=1.$ We find that FEP- is higher than FEP+, implying the transition from the competence state (spiral sink) to the adjacent vegetative region is more likely to occur. Then, when all the noise intensity are fixed at $0.5$ in Fig.6b, with the increase of $\alpha,$ FEP in both directions initially increase, after reaching their maxima, decreases gradually. As shown in Fig.6c, when $\alpha=1,$ in particular with the increase of the noise intensity, FEP- increases, but FEP+ decreases. Then, when all the noise intensity are fixed at $0.5$ in Fig.6d, with the increase of $\alpha,$ FEP in both directions initially increase, after reaching their maxima, drops to [math] gradually. Macroscopically, The behavior of AFEP in both directions are almost same for different scenarios, we omit here.

From the proceeding results on the transition between the vegetative and the competence region, when noise intensity is slightly larger, no matter what escape scenarios we choose, the transition from the high competence region to the low vegetative region is more likely to happen. In the following we explain how to reduce this probability, and how to get the system to reside a longer time in the competence region.

Stochastic basin of attraction

We use now the concept of the basin of attraction of an invariant set [46], i.e., the set whose solution of the system starting from an initial state in it ultimately converges to it. Inspired by refs.29 and 47, we would like to extend the concept of stochastic basin of attraction (SBA), to make a quantification for the stability of a given competence state(spiral sink) in $\mathbb{R}^{2}$.

The stochastic basin of a closed set $K$ is denoted by $\mathfrak{B}(K)\triangleq K\bigcup M$. We employ MFET $u(k,s)$ to quantify the closed set $K$, $K=\{(k,s)\in B|u(k,s)\geq u^{*}\},$ in other words, the stochastic solution starting in $K$ and stays there for a finite time, and $B$ is the selected deterministic basin of attraction of an invariant set. Thus, in the sense of stochastic dynamics, the closed set $K$ can be taken as a stochastic invariant set. It is the collection of all initial states $(k_{0},s_{0})\in\mathbb{R}^{2}$ satisfying: $M=\{(k_{0},s_{0})\in K^{c}|p(k_{0},s_{0})>p^{*}\}$, where $p(k_{0},s_{0})$ represent the probability for all initial points returning to the closed set $K$ from outside of $K$.

Here we suppose that the total noise intensities are fixed, and select the set $B\!=\!(0.115,0.253)\!\times\!(1.06,3.18)$. According to MFET under three fluctuations as listed in Fig.7, we choose $u^{*}=0.0038$ and determine $K\!=\!(0.15295,0.21505)\!\times\!(1.06,3.18).$ The colorful region is the SBA for the competence state (spiral sink). The yellow area stands for the stochastic invariant set $K$ in which the initial points will stay a longer time. Undoubtedly, the competence state (spiral sink) is included in $K.$ The set which has a probability with $p^{*}=0.6$ returns to $K$ is labeled in red and green colors, and the probability with $p^{*}=0.8$ in red color. Fig.7 exhibits a pictorial comprehension for the basin stability of competence state, and there is an expansion in SBA with $\alpha$ decreases. That is, the high competence state is more stable under fluctuations of Brownian motion and $\alpha$-stable Lévy motion with a smaller $\alpha$. In other words, Lévy motion with a larger $\alpha$ is more likely to induce the transition to the low vegetative state.

Discussion

To sum up, we have investigated the impacts of stochastic fluctuations on the evolutionary process of the MeKS network system, where (Gaussian) Brownian and (non-Gaussian) $\alpha$-stable Lévy motion are both incorporated. Two deterministic quantities, the mean first exit time (MFET) and the first escape probability (FEP) from a microscopic perspective, the average first exit time (AMFET) and the average first escape probability (AFEP) from a macroscopic perspective, are employed to analyse the transitions between the vegetative and the competence region in the networks. As we have known that the effects of Lévy and Brownian motions on a one-dimensional gene expression model have recently been investigated [29], we here study a more general two-dimensional MeKS model describing the competence development in the B. subtilis under the influence of Lévy as well as Brownian motion. However, to conduct such a simulation in the two-dimensional model with Lévy motions is such a more challenging but difficult problem. Nevertheless, we have managed to calculate these deterministic quantities by accurate and efficient numerical methods [31, 13, 32, 34, 33, 35].

The noise intensities $\sigma_{k},\sigma_{s},\varepsilon_{k},\varepsilon_{s}$, the Lévy motion index $\alpha$ (relating with the jump magnitude and its frequency) have important and subtle influences on MFET and FEP. We have used the relative contribution factor (RCF) $\lambda$, the ratio between non-Gaussian and Gaussian noise strengths, to implement optimal control in the transition between the two states, when the total quantities of the noise intensities are fixed. We here found that MFET as a function of $\lambda$ shows a maximum, signifying the emergence of the noise-enhanced stability phenomenon for the low vegetative state.

For different initial concentrations, the behavior of the escape probability from the vegetative region to the adjacent competence region is diverse, depending on the influences of the noise intensities and Lévy motion index $\alpha$. The $\alpha$-stable Lévy motion can induce a transition between the vegetative region and a non-adjacent competence region, while Brownian motion can not implement it, due to its continuous paths. Pictorial representations of FEP versus RCF indicate that, for the adjacent scenario, a larger $\alpha$ (a smaller jump magnitude with a higher frequency) makes it happen at the highest probability. While for the non-adjacent scenario, $\alpha_{0}\approx 0.5$ (a larger jump magnitude with a lower frequency) is the best choice. Also, we have surveyed the behavior of escape probability for the transition between the vegetative and the competence region with respect to the two escape scenarios (adjacent, non-adjacent).

Additionally, we have explored the stochastic basin of attraction for the competence state from a geometric viewpoint. Through a schematic presentation of the influence of the Lévy motion index $\alpha$ on the basin stability of the competence region in our stochastic MeKS model, we have detected that the high competence state is more stable under combined fluctuations of Brownian motion and $\alpha$-stable Lévy motion with a smaller $\alpha$.

Furthermore, various diseases pose a tremendous threat to people’s health. Cancer research has become a major challenge in medicine and biology. In complex diseases networks, the gene expression noise induces effects on the choices between the normal and the malignant phenotype of the cell during disease progression. Here, we hope that our exploration will offer a new guidance for the prevention and therapy of complex diseases system with a normal state, a pre-disease state and a disease state.

Methods

The dynamical model of competence induction. To understand how the MeKS network structure determines the dynamics of competence, Süel et al [18] built a mathematical model constrained by experimental observations. This model can be reduced to a system of two ordinary differential equations incorporating both the direct positive and the ComS-regulated negative feedback loops of ComK.

[TABLE]

For convenience, we name it as MeKS model and write $f(k,s)={{a}_{k}}+\tfrac{{{b}_{k}}{{k}^{n}}}{k_{0}^{n}+{{k}^{n}}}-\tfrac{k}{1+k+s}$, and $g(k,s)=\tfrac{{{b}_{s}}}{1+{{(k/{{k}_{1}})}^{p}}}-\tfrac{s}{1+k+s}.$ In the MeKS model, the symbols $k$ and $s$ stand for the concentrations of ComK and ComS protein, respectively. $a_{k}$ and $b_{k}$ denote the basal and fully activated rates of ComK production, and $k_{0}$ the concentration of the ComK protein needed for $50\%$ activation. The cooperativities of ComK auto-activation and ComS repression are prescribed by the Hill coefficients $n$ and $p$, respectively. $b_{s}$ is the maximum expression rate of the ComS protein, and $k=k_{1}$ has the half-maximal rate. The ComS protein competitively inhibits the degradation activity between the MecA-ClpP-ClpC complex and the ComK protein, which leads to the nonlinear degradation terms described at the end of the MeKS model.

It is seen from their nullclines in Fig.1b that there are three equilibria: Nodal sink $(k_{l},s_{l})\approx(0.015262,2.1574)$ (the low stable ComK vegetative state), saddle $(k_{u},s_{u})\approx(0.08568,2.2469)$ (the unstable medium ComK state), and spiral sink $(k_{h},s_{h})\approx(0.15732,1.5781)$ (the high stable ComK competence state).

The stochastic dynamical model of competence induction. To investigate effects of stochastic fluctuations on the transitions between the vegetative and the competence region, our stochastic MeKS model is described as

[TABLE]

Two kinds of noises, Gaussian and non-Gaussian noises, are incorporated in the MeKS model. The symbols $\sigma_{k},\sigma_{s}$ represent the intensities of the Gaussian noise, and $\varepsilon_{k},\varepsilon_{s}$ the intensities of the non-Gaussian noise. The symbols $B_{t}^{1}$ and $B_{t}^{2}$ are two independent Brownian motions, $L_{t}^{1}$ and $L_{t}^{2}$ are two independent $\alpha$-stable Lévy motions, which is a pure jump process and independent of Brownian motions. Due to the vital role of the ComK protein in the competence development, to implement the optimal control of noises for transitions between the vegetative and the competence region, we employ the relative contribution factor $\lambda$, which is the ratio of the intensity of Gaussian noise and non-Gaussian noise, i.e., $\lambda=\frac{\varepsilon_{k}}{\sigma_{k}}$, where $\sigma_{k}$ belongs to $(0,1]$, with the added restrictions $\varepsilon_{k}=\varepsilon_{s},\sigma_{k}=\sigma_{s}$ and $\varepsilon_{k}+\sigma_{k}=1$.

Brief introduction for Lévy motion. A Lévy motion $L_{t}$ is a stochastic process which has stationary and independent increments, and is regarded as an appropriate model for non-Gaussian processes with jumps and heavy tails [48]. In detail, for any time $t_{1},t_{2}$ with $t_{1}<t_{2},$ $L_{t_{2}}-L_{t_{1}}$ and $L_{t_{2}-t_{1}}$ have the same distribution. Namely, the distribution of $L_{t_{2}}-L_{t_{1}}$ only depends on the time difference $t_{2}-t_{1}.$ And for any partition $0=t_{0}<t_{1}<...<t_{n}=t,L_{t_{i}}-L_{t_{i-1}}$ are independent. The sample paths of Lévy motion are continuous from the right and have left limits at every time (càdlàg)[5, 6]. Lévy motion, as a non-Gaussian stochastic process, is a generalization of Gaussian Brownian motion whose sample paths are continuous in time in the common sense.

In our present work, the $\alpha$-stable Lévy motion with a triplet $(0,0,\nu_{\alpha})$, i.e., a pure jump motion is considered. The parameter $\alpha$ is called Lévy motion index, or stability index. Specially, it is a Brownian motion when $\alpha=2.$ For $0<\alpha<2,$ the jump measure in two-dimensional case is defined as $\nu_{\alpha}=\frac{C_{\alpha}dy}{\|y\|^{2+\alpha}},$ where $C_{\alpha}=\frac{\alpha}{2^{1-\alpha}{\pi}}\frac{\Gamma(1+\frac{\alpha}{2})}{\Gamma(1-\frac{\alpha}{2})}$ . The $\alpha$-stable Lévy motion has a larger jump magnitude with a lower jump frequency for $0<\alpha<1$, while it has a smaller jump magnitude with higher jump frequencies for $\alpha$ closes to 2. For more details, we refer readers to an uncomplete list of references [5, 48, 6, 49] etc.

Mean first exit time (MFET) and average mean first exit time (AMFET). The MFET and the AMFET, are utilised to describe the stability for certain stable state of the stochastic MeKS model. First, for a bounded domain $D$ enclosing certain equilibrium, for instance the vegetative state, we define the first exit time as $\tau_{D}(\omega)=\inf\{t\geq 0,X_{0}=x,\leavevmode\nobreak\ X_{t}\in D^{c}\},$ where $D^{c}$ denotes the complement set of $D$ in $\mathbb{R}^{2}$. Then the mean exit time is denoted as $u(x)\triangleq\mathbb{E}\tau_{D}(\omega)\geq 0$. It is the mean residence time in certain domain $D$ before exiting to another regime and used to our analysis from microscopic perspective. According to ref.6, the MFET $u(k,s)$ in our model can be solved by the following differential-integral equations with the Dirichlet exterior boundary condition

[TABLE]

where $D=(a,b)\times(c,d)$ with real numbers $a,b,c,d$. Here, due to the independent properties of the Brownian and Lévy motion, the jump measure in two-dimensional case is divided into the summation of two one-dimensional jump measure [50]. Then the generator $\mathscr{A}$ is

[TABLE]

Secondly, the AMFET on $D$ is defined as $\bar{u}(D)=\frac{1}{|D|}\int\int_{D}u(k,s)dkds$ for our analysis from the macroscopic perspective [6, 35].

First escape probability (FEP) and average first escape probability (AFEP). The escape probability here are employed to characterize the likelihood of the transition between the low vegetative and the high competence region for the MeKS system. The FEP $p(k,s)$ can be solved by the following differential-integral equations with the Balayage-Dirichlet exterior boundary condition

[TABLE]

Here $\mathscr{A}$ is the same generator as (5).

The AFEP on $D$ is defined as $\bar{p}(D)=\frac{1}{|D|}\int\int_{D}p(k,s)dkds$ [6, 35]. We use the FEP to examine the transitions from a microscopic viewpoint, and the AFEP from the macroscopic perspective.

It has recently been examined that the effects of Lévy and Brownian motions on a one-dimensional gene expression model [29]. Here, we study a more general stochastic two-dimensional MeKS model describing the competence development in the B. subtilis. However, it turns into a more challenging but difficult problem to conduct such a simulation in the two-dimensional model with Lévy motions. Nevertheless, we have managed to apply recent results to calculate these deterministic quantities by accurate and efficient numerical methods [31, 13, 32, 34, 33, 35]. Fig.2 exhibits a visualized representation for the MFET and the FEP in the two-dimensional stochastic MeKS model (3).

Numerical methods. To solve the MFET and the FEP efficiently and accurately, Gao et al [31] constructed an effective numerical method to solve the aforementioned differential-integral equations. They used the central difference for derivatives and a modified trapezoidal rule to discretize the differential-integral equations [31]. Their scheme can be modified for our model with two transformations, i.e., $k=\frac{b-a}{2}z+\frac{b+a}{2}$ for $z\in(-1,1)$, and $s=\frac{d-c}{2}w+\frac{d+c}{2}$ for $w\in(-1,1)$, correspondingly, $k^{\prime}=\frac{b-a}{2}\hat{k}$ and $s^{\prime}=\frac{d-c}{2}\hat{s}$. We let $v(z,w)=u(\frac{b-a}{2}z+\frac{b+a}{2},\frac{d-c}{2}w+\frac{d+c}{2})$, and address the integral term in equation (5) by the Cauchy principal value integral. Then the generator $\mathscr{A}$ can be rewritten as

[TABLE]

Here $\phi(z,w)$ stands for the right-hand side term occurring at the differential-integral equation. For the equation (4), $\phi(z,w)=-1.$ While for the equation (6), the expression of $\phi(z,w)$ varies with the different escape scenarios. In detail, first, for the adjacent scenario, if the system escapes from region $D=(a,b)\times(c,d)$ to an adjacent region $E=(b,i)\times(c,d),$ then $\phi(z,w)=\frac{\varepsilon_{k}^{\alpha}C_{\alpha}}{\alpha}\{(i-\frac{b-a}{2}z-\frac{b+a}{2})^{-\alpha}-[\frac{b-a}{2}(1-z)]^{-\alpha}\},$ especially, for $i=\infty,$ $\phi(z,w)=-\frac{\varepsilon_{k}^{\alpha}C_{\alpha}}{\alpha}(\frac{2}{b-a})^{\alpha}\frac{1}{(1-z)^{\alpha}}.$ Reversely, if the system escapes from region $D=(a,b)\times(c,d)$ to an adjacent region $E=(h,a)\times(c,d),$ then $\phi(z,w)=\frac{\varepsilon_{k}^{\alpha}C_{\alpha}}{\alpha}\{(\frac{b-a}{2}z+\frac{b+a}{2}-h)^{-\alpha}-[\frac{b-a}{2}(1+z)]^{-\alpha}\}.$ Secondly, for the non-adjacent scenario, if the system escapes from region $D=(a,b)\times(c,d)$ to a non-adjacent region $E=(h,i)\times(c,d),$ then $\phi(z,w)=\frac{\varepsilon_{k}^{\alpha}C_{\alpha}}{\alpha}\{(i-\frac{b-a}{2}z-\frac{b+a}{2})^{-\alpha}-(h-\frac{b-a}{2}z-\frac{b+a}{2})^{-\alpha}\}.$ Reversely, if the system escapes from region $D=(a,b)\times(c,d)$ to a non-adjacent region $E=(h,i)\times(c,d),$ then $\phi(z,w)=\frac{\varepsilon_{k}^{\alpha}C_{\alpha}}{\alpha}\{(\frac{b-a}{2}z-\frac{b+a}{2}-h)^{-\alpha}-(\frac{b-a}{2}z-\frac{b+a}{2}-i)^{-\alpha}\}.$

Acknowledgements

We would like to thank Dongfang Li, Xiao Wang and Hansen Ha for discussions about computation. We are grateful to Aining Bai in Institute of Botany (CAS), Liang Wang, Mengli Hao, Jian Ren, Tao Jiang, Xinyong Zhang, Yanjie Zhang, Ziying He for helpful discussions. This work was partly supported by the NSFC grants 11531006, 11371367 and 11271290.

Author contributions statement

F.W. and J.D. conceived the research and wrote the first draft of the manuscript. F.W. performed computations. X.C. and X.L. conducted and verified the program designs. Y.Z. and J.K. analysed the results and concepts development. All authors conducted research discussions and reviewed the manuscript.

Additional information

Competing financial interests: The authors declare no competing financial interests.