The gene regulatory network ( Figure 1 ) of IGC contains 17 genes. Through the collection of our simulation trajectories and associated statistics, we can quantify the landscape on a 17-dimensional probability distribution. The related potential landscape U can be defined as U = -lnP_ss (41, 42). P_ss is the probability of the steady-states, with’ss’ being an abbreviation of steady state. It is difficult to visually display a 17-dimensional landscape, so we take two genes or dimensions (HIF - 1α and COX2, two genes very crucial to IGC) to visualize the landscape clearly. In Figure 2 , the X-axis shows the expression level of hypoxia-inducible factor-1α (HIF-1α), which plays the role of an ‘angiogenic switch’ in the hypoxia microenvironment in many types of tumors, including IGC (43). The Y-axis shows the expression level of Cyclooxygenase-2 (COX2), which plays a crucial role in cancer development and clinical metastasis (44). We can also choose two other genes (such as EGFR and VEGF), as shown in the Supplementary Material . Here, there are still three stable states. There are some changes in the landscape topography, as different genes have different values for the three stable points.

On the IGC landscape, there are three stable state attractors, which are the normal, AG, and cancer (IGC) states, respectively. The definition of these three stable states is based on the biological functions and gene expression levels of the 17 genes in the gene regulatory network. The parameters we set in our model depend on the simulated gene expression levels of the 17 genes, and are all in agreement with the trends seen in the clinical data at different cancer stages (different states). P53 and APC are tumor suppressor genes that have high expression levels in the normal state but low expression levels in cancer cells. The other 15 genes have high expression levels in cancer and low expression levels in normal cells. The simulated gene expression levels of the 17 genes are all consistent with the experimental literature (details can be seen in Supplementary Material Table S3 ).

In Figure 2 , the genes HIF-1α and COX2 have high expression levels in the cancer state, low expression levels in the normal state, and intermediate expression levels in the intermediate state (AG state) (45). The barrier height between the normal and AG states is relatively low, which indicates that AG infection and recovery is relatively achievable. The barrier height between AG and gastric cancer states is much higher, which indicates that gastric cancer formation and recovery is much more difficult. The barrier height can quantify the difficulty of attractor transfer from one state to another. The barrier height from gastric cancer to the AG state is very high, which shows why gastric cancer is so difficult to cure (or reverse).

To describe the gastric cancer progression process quantitatively, we quantified the dominant paths from normal to AG to gastric cancer states using previously-explored approaches (path integral approaches) (46, 47). In Figure 2B , we can see the dominant paths colored in red, blue, yellow, and purple, respectively, as the dominant paths from the normal to the AG state, from the AG to the normal state, from the AG to the gastric cancer state, and from the gastric cancer to the AG state, respectively. These dominant paths are separated and irreversible (48) as the rotational flux force (as part of the driving force) in addition to the gradient force makes the dominant path separate from the gradient direction of the potential. In our model, the driving force can be mathematically decomposed into two directions, the rotational flux force and the gradient force of the potential landscape. The green arrows represent the rotational flux force and the white arrows represent the gradient force direction. The dominant paths through the normal state to the AG and gastric cancer states are irreversible, which can help us understand why the processes of IGC formation and IGC treatment are separate and irreversible biological processes.

To investigate the influence of H. pylori on IGC, we performed simulations to observe the cancer progression on the landscape. We used a term in the ODEs to simulate different degrees of H. pylori infection. The term F(x_i ) can be rewritten as F′(x_i ) = F(x_i ) + H_i (i = 1,2…,17). The term H_i is used to denote the degree of Helicobacter pylori infection on the related gene expression level. The value of H is set according to the experiments. If the gene expression level is increased from the H. pylori infection, the value of H will be >0. In the opposite case, H < 0.

For Figure 3 , we chose genes hTERT and MYC to show the landscape layers with variations associated with the development of IGC under the effects of H. pylori infection. H. pylori infection can result in the gene expression levels of both genes being increased. From Figure 3 , we can see that when the term H = 0, the normal, AG, and gastric cancer states are visible on the first layer of the landscape. A value of H = 0.05 indicates infection. There are normal, AG, IM, and gastric cancer states on the second layer of the landscape. The probability of the AG state is dominant, indicating that the H. pylori infection accelerates the development of AG. When the term H = 0.1, there are normal, IM, and gastric cancer states on the third layer of the landscape. The AG state has disappeared and the IM state is dominant, which indicates that the H. pylori infection worsens gastritis and causes further change into intestinal metaplasia. When the term H = 0.2, there are normal and gastric cancer states on the fourth layer of the landscape. The IM and gastric cancer states gradually converge and merge into one. The cancer state becomes the dominant state. When the term H = 0.6, there is only one gastric cancer state on the fifth layer of the landscape. The normal state disappears and the cancer state is dominant. It is impossible for a patient to recover to their normal state while suffering from H. pylori infection can lead to the aggravation of AG, and then the appearance of the IM state. The IM state can be considered very close to the cancerous state during IGC development. When the H. pylori infection becomes more and more serious, the AG state disappears and the IM state becomes dominant, finally leaving only one gastric cancer state.

This simulation is on the epigenetic level to illustrate the progression and development of IGC when one gets infected with H. pylori. The regulations of the network do not change with this series of variations in the landscape. The effect of H. pylori infection results in variations in the IGC landscape. Depending on the degree of H. pylori infection, which is becoming more and more serious, landscape development is moving on the cancer direction. As AG and IM are seen in the development of the landscape, which demonstrates the dynamic process of IGC development as the degree of H. pylori infection changes.

To further investigate the key regulations crucial to IGC therapy or prevention, we apply a global sensitivity analysis method to the landscape model. Each gene regulation or protein concentration can contribute to system dynamics. A small change in the regulatory strength of one gene in the gene regulatory network can lead to the whole landscape topography varying accordingly. As the barrier heights between the biological states can quantify the difficulty associated with transferring between the states, we calculate the variation in the barrier height when regulation strength is changed. The greater the variation, the more sensitive the regulation is.

Figure 4 shows the global sensitivity analysis for IGC. We selected the top 10 most sensitive regulations, which are shown in Figure 4A . When we changed the regulation strength to 0.9 of the original regulation strength, these regulations showed the most significant variations in ΔU_ng and ΔU_gn . ΔU_ng is the variation of the barrier height from the normal state to the AG state. ΔU_gn is the variation of the barrier height from the AG state to the normal state. From Figure 4A , we can see that the values of ΔU_ng changed most significantly are for regulating RAS → HIF-1α and ZEB → TGF-β. The higher value of ΔU_ng indicates that it is more difficult for the cells to transform from the normal state to the AG state than before. This is because the barrier height between the normal and the AG state is much higher than before. The variations of ΔU_ng are most significant when the regulation strengths of the RAS → HIF-1α and ZEB → TGF-β are varied. This type of variation can be used in gastritis prevention as the cell transformation to the AG state becomes more difficult. When the regulation strength is reduced to 0.9 of the original value, the activation of the genes HIF-1α and TGF-β decreases and the concentrations of HIF-1α and TGF-β decrease accordingly. Experiments show that the expression levels of HIF-1α in gastric cancer patients are higher than those in healthy subjects (25). HIF-1α participates in the activation of numerous target genes to adapt to the hypoxic environment (49), which leads to gastric cancer development. We changed the regulation strength, which can cause the concentration of HIF-1α to decrease. This will inhibit the transcription of those target genes and reduce the ability of the cells to adapt to the hypoxic environment, which leads to higher ΔU_ng . Therefore, reducing the regulation strength to 0.9 times the original regulation strength can inhibit the cells from adapting to the hypoxic environment, which can prevent the gastritis cells from developing into gastric cancer. Studies show that inhibiting TGF-β expression can serve as a potential therapeutic target or a biomarker for gastric cancer treatment (50, 51). Therefore, the decreased regulation strength can reduce the transcription of TGF-β, which increases in ΔU_ng and inhibits gastritis from developing into gastric cancer.

Figure 4B shows the top 10 most sensitive regulations when we changed the regulation strengths to 0.9 times their original values. These regulations show the most significant variations of ΔU_gc and ΔU_cg . A higher value of ΔU_gc compared to before indicates that switching from the AG state to the gastric cancer state becomes more difficult, because the barrier height between the AG and gastric cancer states is much higher. Reducing the regulation strength can inhibit the expression level of RAS, which increases in ΔU_gc and helps inhibit gastric cancer progression and development. Such variation can be used in gastric cancer prevention, as the transformation of the cell from the AG state to the gastric cancer state becomes more difficult. When the regulation strength is reduced to 0.9 times the original value, the activation of the RAS gene decreases and the concentration of RAS also decreases accordingly. There are many studies showing that K-ras is associated with the development and progression of gastric cancer. Overexpression of K-ras can increase the risk of gastric cancer development (52) (K-ras is a protein of the RAS family). Reducing the regulation strength can inhibit the expression level of RAS, which increases in ΔU_gc and helps inhibit gastric cancer progression and development.

Figure 4 displays the most sensitive regulations on global topography in terms of the barrier height between normal and cancer states. When the regulation strengths of P53 → RAS, c-myc→ RAS, HIF-1α → RAS and TGF-β → RAS are reduced to 0.9 of the original values, the expression of RAS will decrease accordingly. Studies show that K-ras regulates cell survival, motility, proliferation, angiogenesis, and metastasis (53). Therefore, this participates in gastritis and gastric cancer formation and metastasis. When the regulation strengths of RAS → HIF-1α, VEGF → HIF-1α, and Bcl-2 → HIF-1α are reduced to 0.9 of the original values, the expression of HIF-1α is reduced accordingly. The gene HIF-1α can activate the transcription of many target genes to adapt to the hypoxic environment of cancer cells (54). The overexpression of the gene HIF-1α can induce cancer cell development. When the regulation strengths ofP53 → TGF-β and ZEB → TGF-β are reduced to 0.9 of the original values, TGF-β expression decreases. TGF-β promotes cancer-related characteristics in most gastric cancer cell lines (55). In TGF-β and HIF-1α gene expressions appear twice as that of the target genes in Figure 4A . Ras expression appears thrice in Figure 4 , and the expressions of HIF-1α, P53, and hTERT genes appear twice as that of the target genes. We can pay more attention to these genes in designing strategies in clinical experiments or trials to prevent gastritis or gastric cancer formation.

In Figure 3 , we can see that when a patient is infected with H. pylori (term H = 0.05), four stable states emerge, which are the normal, AG, IM, and gastric cancer states. To figure out which regulations are more sensitive to IGC with H. pylori infection, we performed a global sensitivity analysis on this condition.

Figure 5 shows global sensitivity analysis for IGC with H. pylori infection. We reduced the regulation strength to 0.9 of the original value. Figure 5A displays the top 10 regulations most sensitive to variations in ΔU_gm and ΔU_mg . ΔU_gm is the variation in the barrier height from the AG state to the IM state. ΔU_mg is the variation in the barrier height from the IM state to the AG state. The regulations for HIF-1α→c-myc and CDK2→c-myc cause the most significant changes in ΔU_gm . The value of ΔU_gm becoming higher indicates that it becomes more difficult for cells to transform from the AG state to the IM state as the barrier height is higher. This type of variation can be used to prevent cell transformation from the AG state to the IM state. As the activation strengths of HIF-1α→c-myc and CDK2→c-myc are decreased, the expression of c-myc will be reduced accordingly. Gene c-myc has been studied as a biomarker with which to identify H. pylori infection (56). Gastric cancer treatment and gastric cancer progression are complicated by aberrant expressions of c-myc (57). Therefore, inhibiting the expression of c-myc will benefit the treatment of IGC with H. pylori infection and prevent transformations from AG to an IM state.

Figure 5B displays the top 10 regulations most sensitive to variations in ΔU_mc and ΔU_cm . ΔU_mc is the variation of the barrier height from the IM state to the gastric cancer state. ΔU_cm is the variation of the barrier height from the gastric cancer state to the IM state. The regulation RAS→ TGF-β caused the most significant changes in ΔU_mc . The value of ΔU_mc becoming higher indicates that it becomes more difficult for cells to transform from the IM state to the gastric cancer state because the barrier height between the two states is higher. As the activation strength of RAS→ TGF-β is decreased, the expression of TGF-β will be reduced accordingly. TGF-β can trigger epithelial–mesenchymal transition (EMT) markers, which are crucial for canceration and metastasis. H. pylori infection can induce TGF- β to trigger the EMT pathway. When H. pylori is eradicated, TGF- β is inhibited from triggering the EMT pathway (58). TGF- β is a key gene in gastric cancer prevention and treatment, which we have identified in this study. C-myc and TGF- β are vital for the treatment and prevention of IGC from H. pylori infection.

From Figure 4 and Figure 5 we can see that the key regulations are different depending on whether the H. pylori is infected or uninfected. When infected with H. pylori, another state (IM state) emerged on the landscape, which is different from the IGC landscape without H. pylori infection. The c-myc gene is essential as it appears in the two key regulations (HIF-1α→c-myc and CDK2→c-myc). The c-myc gene is a biomarker to identify H. pylori infection in clinical trials (56). We can take other regulations such as TGF-β→ IL-1β and TGF-β→ ZEB, which are more sensitive in clinical experiments. The gene TGF-β is vital for both H. pylori infected and uninfected as it is sensitive to the two conditions when the cell states switch from AG (or IM) state to cancer state. TGF-β plays a critical role in cancer metastasis (58). TGF-β appears thrice in the top 10 regulations in Figure 5B , while TNF-α and c-myc genes appear twice. We should take these genes into consideration in designing strategies in clinical experiments for preventing of gastric cancer with H. pylori infection.

A gene regulatory network consists of genes and gene regulatory relationships, represented by nodes and edges in the network, respectively. We use an n-component vector x = (x ₁, x ₂,…, x_n ) to quantify the genes in the network. Here n is the number of genes in the network and x_i (i =1,2,…, n) denote the expression levels (or protein concentrations) of the corresponding genes. A system of ordinary differential equations, written in the compact form x ˙ = F ( x ) , can be employed to study the deterministic dynamics of the network, where F ( x ) denotes the driving force of the deterministic dynamics.

In biological systems, stochastic fluctuations of internal or external origins are ubiquitous and may have a significant impact on the dynamics of the system. To incorporate the effects of stochastic fluctuations, a stochastic force ξ(t) may be attached to the ordinary differential equations x ˙ = F ( x ) . This leads to a stochastic differential equation of the form x ˙ = F ( x ) + ξ ( t ) , also known as the Langevin equation. The stochastic force ξ(t) modeling random fluctuations is assumed to be Gaussian white noise in time, with the mean〈 ξ(t) 〉=0 and the correlation 〈 ξ(t)ξ ^T (t') 〉=2Dδ(t−t'). Here D is the diffusion matrix characterizing the fluctuation strength.

An equivalent description of the Langevin dynamics is in terms of the probability distribution P( x , t), whose time evolution is governed by the corresponding Fokker-Planck equation: ∂P/∂t = -∇ · [ F P - ∇· ( D P)]. It can also be written as ∂P/∂t = - ∇ · J , with J denoting the probability flux. The steady state characterized by ∂P_ss /∂t = ∇ · J _ss = 0 is of particular interest. In an equilibrium system, the probability flux at the steady state vanishes, i.e. J _ss = 0. In a non-equilibrium system, there is in general a nonvanishing probability flux at the steady state, i.e. J _ss  ≠ 0, which signifies the time-irreversible nature of the non-equilibrium steady state. From the expression J _ss = F P_ss – D · ∇P_ss , the driving force F can be written in the landscape-flux decomposition form (59): F = - D · ∇U + J _ss /P_ss . Here – D · ∇U is the part of the driving force contributed by the gradient of the potential landscape U = -ln P_ss , and J _ss /P_ss is the other part contributed by the probability flux that is associated with the nonequilibrium nature of the system.

Self-consistent mean field (60) serves as an effective approximation method of solving Fokker-Planck equations with a large number of variables. In this approximation, the joint probability distribution of all the variables is substituted by the product of marginal probability distributions of each variable, namely, P(x ₁, x ₂,…, x_n , t) ~ Π _i P (x_i , t), so that the latter can be solved in a self-consistent manner. The dimensionality of the problem is reduced significantly from mⁿ to m × n, where m is the number of possible values each variable may take. This makes the computations much more feasible.

A further approximation is invoked to simplify the problem, which postulates P(x ₁,x ₂,…,x_n , t) with the form of a multivariate Gaussian distribution. When the magnitude of the diffusion matrix D is small, the equations governing the mean vector x ¯ ( t ) and the covariance matrix σ (t) of the Gaussian distribution are given by:

Here the matrix A has elements A i j ( t ) = ∂ F i ( x ¯ ( t ) ) ∂ x ¯ j ( t ) . Given the self-consistent mean field approximation, only the diagonal elements of σ(t) need to be considered. The combination of the self-consistent mean field approach and the Gaussian distribution approximation leads to the following form of probability distribution evolution for each x_i :

The steady-state probability distribution of a monostable system with one fixed point can be easily obtained on the basis of Eq.(5) as a single Gaussian distribution. In a multistable system with more than one fixed point, the steady-state distribution may be constructed as a combination of multiple Gaussian distributions with the form P_ss ( x ) = Σ_kω_kP_k ( x ), where k labels different fixed points, ω_k represents the weight of each fixed point, and P_k ( x ) is the Gaussian distribution corresponding to each fixed point.

Based on the Onsager-Machlup functional, the transition probability of the Fokker-Planck equation has the following path-integral formulation (46):

In the above, x ₀ denotes the initial state at time t ₀ and x _f represents final state at time t_f . P( x _f , t_f ; x ₀, t ₀) is the transition probability from the initial state to the final state. The notation ∫ D [ x (t)] represents an integral over the all possible paths starting from the initial state x ₀ at time t ₀ and ending at the final state x _f at time t_f . L( x (t)) is the Lagrangian with the expression L ( x ( t ) ) = 1 4 ( x ˙ − F ( x ) ) · D − 1 · ( x ˙ − F ( x ) ) + 1 2 ∇ · F ( x ) . Its time integration gives the action S[ x (t)] = ∫ L( x (t))dt associated with each path as in classical mechanics. The action S[ x (t)] determines the probability weight e ^{- S [ x ( t )]} contributed by the corresponding path. The summation (or integration) of these probability weights over all the paths gives the transition probability. Since the contribution of each path has the exponential form e ^{- S [ x ( t )]}, the dominant path with maximum probability is the path with minimum action, which can be determined by the variational principle δS[ x (t)] = 0 and the resulting Euler-Lagrange equation. For non-equilibrium systems the existence of nonvanishing probability flux J _ss cannot be ignored. As a consequence, the dominant kinetic paths in non-equilibrium systems are separated and irreversible (46).