Augmented Variational Superposed Gaussian Approximation for Langevin Equations with Rational Polynomial Functions

TL;DR

This paper introduces an augmented variational superposed Gaussian approximation method that efficiently solves Langevin equations with rational polynomial drift functions, reducing computational costs compared to Monte Carlo simulations.

Contribution

The paper proposes an augmented VSGA method using Padé approximant to handle non-polynomial drift functions in Langevin equations, improving efficiency and accuracy.

Findings

Accurately solves Langevin equations with rational polynomial drift functions.

Reduces computational cost compared to Monte Carlo methods.

Successfully applied to chaotic, genetic, and bistable systems.

Abstract

Reliable methods for obtaining time-dependent solutions of Langevin equations are in high demand in the field of non-equilibrium theory. In this paper, we present a new method based on variational superposed Gaussian approximation (VSGA) and Pad\'e approximant. The VSGA obtains time-dependent probability density functions as a superposition of multiple Gaussian distributions. However, a limitation of the VSGA is that the expectation of the drift term with respect to the Gaussian distribution should be calculated analytically, which is typically satisfied when the drift term is a polynomial function. When this condition is not met, the VSGA must rely on the numerical integration of the expectation at each step, resulting in huge computational cost. We propose an augmented VSGA (A-VSGA) method that effectively overcomes the limitation of the VSGA by approximating non-linear functions with the Pad\'e approximant. We apply the A-VSGA to two systems driven by chaotic input signals, a stochastic genetic regulatory system and a soft bistable system, whose drift terms are a rational polynomial function and a hyperbolic tangent function, respectively. The numerical calculations show that the proposed method can provide accurate results with less temporal cost than that required for Monte Carlo simulation.

Figures

Equations

$$\frac{d{x}_i}{dt}=f_i(\mathbf{x},t)+\sum _{j=1}^{N_g}{g}_{ij}(\mathbf{x},t){\xi }_j(t),(i=1,2,...,N),$$

$$\frac{\partial }{\partial t}P(\mathbf{x};t)=\hat{L}(\mathbf{x};t)P(\mathbf{x};t),$$

$$\hat{L}(\mathbf{x};t)=-\sum _i\frac{\partial }{\partial x_i}\left[f_j(\mathbf{x})+\sum _{k,j}{g}_{kj}(\mathbf{x};t)\frac{\partial }{\partial x_k}{g}_{ij}(\mathbf{x};t)\right]+\sum _{i,j}\frac{{\partial }^2}{\partial x_i\partial x_j}\sum _k{g}_{ik}(\mathbf{x};t)g_{jk}(\mathbf{x};t).$$

$$R[\Theta ]={\int }_{-\infty }^{\infty }{\left[\hat{L}(\mathbf{x};t)P(\mathbf{x};t)-\Theta (\mathbf{x};t)\right]}^{2}d\mathbf{x}.$$

$$\tilde{R}[\Theta ]={\int }_{-\infty }^{\infty }{\left[\hat{L}(\mathbf{x};t)P(\mathbf{x};t)-\Theta (\mathbf{x};t)\right]}^{2}d\mathbf{x}+\lambda (t){\int }_{-\infty }^{\infty }\Theta (\mathbf{x};t)d\mathbf{x}.$$

$${\int }_{-\infty }^{\infty }\delta \Theta \left[\hat{L}(\mathbf{x};t)P(\mathbf{x};t)-\Theta (\mathbf{x};t)+\lambda (t)\right]d\mathbf{x}=0.$$

$$P(\mathbf{x};t)=P(\mathbf{x};\mathbf{\theta }(t))=\sum _{m=1}^{N_B}{r}_m(t)\exp \left[-{\mathbf{x}}^{\top }{\mathbf{A}}_m(t)\mathbf{x}+{\mathbf{b}}_m^{\top }(t)\mathbf{x}\right],$$

$$K=\frac{N_B(N+1)(N+2)}{2}.$$

$${\int }_{-\infty }^{\infty }\frac{\partial P(\mathbf{x};\mathbf{\theta })}{\partial {\theta }_l}\left[\hat{L}(\mathbf{x};t)P(\mathbf{x};\mathbf{\theta })-\Theta (\mathbf{x};\mathbf{\theta },\dot{\mathbf{\theta }})+\lambda (t)\right]d\mathbf{x}=0,(l=1,2,...,K).$$

$$f(x)=\frac{g(x)}{h(x)},$$

$$\frac{\partial }{\partial t}P(x;t)=\Theta (x;t)=-\frac{\partial }{\partial x}\frac{g(x)}{h(x)}P(x;t)+D\frac{{\partial }^2}{\partial x^2}P(x;t).$$

$$0=-h^2(x)\Theta (x;t)-\left[\frac{\partial g(x)P(x;t)}{\partial x}h(x)-\frac{\partial h(x)}{\partial x}g(x)P(x;t)\right]+D{h}^2(x)\frac{{\partial }^2}{\partial x^2}P(x;t).$$

$$R[\Theta ]={\int }_{-\infty }^{\infty }{\left[-\frac{\partial g(x)P(x;t)}{\partial x}h(x)+\frac{\partial h(x)}{\partial x}g(x)P(x;t)+D{h}^2(x)\frac{{\partial }^2}{\partial x^2}P(x;t)-h^2(x)\Theta (x;t)\right]}^{2}dx.$$

$${\int }_{-\infty }^{\infty }\frac{\partial P(x;\theta )}{\partial {\theta }_l}[\hat{L}(x;t)P(x;\theta )-h^4(x)\Theta (x;\theta ,\dot{\theta })+\lambda (t)]dx=0,(l=1,2,\dots ,3N_B),$$

$$\hat{L}(x;t)P(x;\theta )=\frac{\partial g(x)P(x;\theta )}{\partial x}{h}^3(x)-\frac{\partial h(x)}{\partial x}g(x)P(x;\theta )h^2(x)+D{h}^4(x)\frac{{\partial }^2}{\partial x^2}P(x,\theta )).$$

$$V(x)\approx \frac{a_0+a_{1}x+a_2{x}^2+...+a_m{x}^m}{1+b_{1}x+b_2{x}^2+..+b_n{x}^n}=\frac{g(x)}{h(x)},(m\ge 0,n\le 1,m,n\text{are integers}).$$

$$\frac{dx}{dt}=\frac{k_f{x}^2}{x^2+K_d}-k_{d}x+R_{\mathrm{bas}},$$

$$\frac{dx}{dt}=\frac{k_f{x}^2}{x^2+K_d}-k_{d}x+I(t)+R_{\mathrm{bas}}+\sqrt{D}\xi (t),$$

$$\hat{L}(x;t)=-\frac{\partial }{\partial x}\left[\frac{k_f{x}^2}{x^2+K_d}-k_{d}x+I(t)+R_{\mathrm{bas}}\right]+D\frac{{\partial }^2}{\partial x^2}.$$

$$P_{st}(x)=\frac{1}{Z(D)}\exp \left[-\frac{U(x)}{D}\right],$$

$$\frac{d{s}_1}{dt}$$

$$\frac{d{s}_2}{dt}$$

$$\frac{d{s}_3}{dt}$$

$$I(t)=\alpha s_1(\omega t),$$

$$\frac{dx}{dt}=-\frac{x}{2}+\tanh x+I(t)+\sqrt{D}\xi (t),$$

$$\frac{g(x)}{h(x)}=\frac{x+\frac{x^3}{15}}{1+\frac{2x^2}{5}}.$$

$$\hat{L}(x;t)=-\frac{\partial }{\partial x}\left[\frac{x+\frac{x^3}{15}}{1+\frac{2x^2}{5}}-\frac{x}{2}+I(t)\right]+D\frac{{\partial }^2}{\partial x^2}.$$

Full text

I Introduction

In recent years, non-equilibrium systems have attracted much interest by virtue of advancements in non-equilibrium thermodynamics Ritort (2008); Seifert (2012). To investigate the dynamics of non-equilibrium systems, studies have focused on the properties and analytical solutions of Langevin equations Hida (1980); Van Kampen (1992); Coffey and Kalmykov (2004); Fa (2011). For stationary systems, the solutions of Langevin equations do not depend on time and can be calculated analytically for one-dimensional systems. However, systems are driven by external forces in dynamic cases. Except for some particular cases, time-dependent solutions cannot be obtained analytically for such driven systems, not even for one-dimensional cases Risken (1996). Therefore, the time-dependent solutions of Langevin equations with external forces or input signals remain a difficult problem.

Our goal is to propose a reliable numerical method that can be used to solve time-dependent solutions of Langevin equations with driving forces to meet the needs of non-equilibrium theory Ritort (2008); Seifert (2012); Weeks and Gilmer (2007). The moment method is a widely applied approach which computes the time evolution of the mean and the covariance of Langevin equations Rodriguez and Tuckwell (1996); Tuckwell and Jost (2009). Although it is simple and computationally efficient, its accuracy is not satisfactory even for simple bistable models. Direct numerical approaches such as a finite-element method Harrison (1988); Kumar and Narayanan (2006) and a finite-difference method Whitney (1970) were proposed to obtain time-dependent probability density functions (PDFs), but the computational cost associated with the methods was very high. Considering the success of the multiple Gaussian method with the variational principle for time-dependent wave functions in quantum mechanics Heller (1976), several studies have tried to approximate solutions of Langevin equations by a superposition of Gaussian distributions Er (1998); Pradlwarter (2001); Terejanu et al. (2008); Hasegawa (2015). Er Er (1998) obtained PDFs as a superposition of multiple Gaussian distributions with the weighted residual method, but the approach is limited to stationary cases. Pradlwarter Pradlwarter (2001) represented time-dependent solutions by a sum of small elements of Gaussian distributions. A disadvantage of this method is that it must handle the variance and the number of Gaussian distributions. Terejanu et al. Terejanu et al. (2008) proposed a superposition of Gaussian distributions as time-dependent PDFs. Their approach first assumed that the weights for each of the Gaussian distributions are fixed. After calculating the mean and variance of the Gaussian distributions, the weights are optimized by minimizing the squared error. The variational superposed Gaussian approximation (VSGA) is shown to be effective and accurate in the calculation of time-dependent solutions in strongly non-linear systems Hasegawa (2015). When compared with Er (1998); Terejanu et al. (2008); Pradlwarter (2001), the VSGA directly obtains the time evolution of the mean, variance, and weight and does not require extra steps. Still, a drawback of the VSGA is that the expectation of the drift term with respect to the Gaussian distribution should be obtained analytically (e.g., when the drift term is a polynomial function). Otherwise, the expectation would be calculated by time-consuming numerical integration. This requirement greatly restricts the application of the VSGA.

In this paper, we propose an augmented VSGA (A-VSGA) method that can handle systems with rational polynomial functions. Rational functions can be used to represent many dynamical systems such as biochemical reactions, gene regulation, and dose-response. Furthermore, through Padé approximant, non-linear systems with drift terms (irrational functions) can be approximated by rational functions. The use of Padé approximant enables applications of the A-VSGA to such non-linear systems. In many cases, the Padé approximant provides closer approximation than the Taylor series expansion, and it works even when the Taylor expansion diverges at poles. To test the accuracy and effectiveness of the A-VSGA, we first apply the A-VSGA to a stochastic gene regulatory system where the drift term is represented by a Hill equation. The Hill equation is a fundamental equation to describing the cooperative effects of multimers. It has been widely used in systems biology Weiss (1997); Tian (2010); Smolen et al. (2000); Liu and Jia (2004). Next, to investigate its reliability, we apply the A-VSGA with the Padé approximant to a soft bistable system subject to white Gaussian noise. The soft bistable model reflects the dynamics of the membrane potential of a neuron cell and is represented by a hyperbolic tangent function Longtin et al. (1994). As noted above, its drift term can be approximated by the rational polynomial functions via the Padé approximant. Bistability describes the switching dynamics of systems and often occurs in non-linear systems far from equilibrium Wilhelm (2009). Due to non-linearity, many researches have studied bistable systems by simulation only. We demonstrate our method as an effective and accurate approach for solving time-dependent solutions of Langevin equations with a non-polynomial expression of the drift terms.

The paper is organized as follows. In Sec. II, we briefly introduce the VSGA and its shortcomings, and propose an augmented method, the A-VSGA. Then, we investigate the accuracy and effectiveness of the A-VSGA in Sec. III by applying it to a gene regulatory system and a soft bistable system driven by time-dependent input signals. Finally, we discuss the problems in the experiments and summarize our findings in Sec. IV.

II Methods

An $N$-dimensional Langevin equation,

[TABLE]

describes the temporal evolution of $N$ state variables $\bm{x}=(x_{1},x_{2},...,x_{N})^{\top}$. In Eq. (1), $f_{i}(\bm{x},t)$ and $g_{ij}(\bm{x},t)$ denote drift and multiplicative terms, respectively. $\xi_{i}(t)$ is white Gaussian noise with zero mean $\left\langle\xi_{i}(t)\right\rangle=0$, and the correlation is $\left\langle\xi_{i}(t)\xi_{j}(t^{\prime})\right\rangle=2\delta_{ij}\delta(t-t^{\prime})$, where $\delta$ is the Dirac delta function. We interpret Eq. (1) in the Stratonovich sense. At a given time $t$, the state of system is not deterministic due to the noise, $\bm{\xi}(t)$. By repeatedly solving Eq. (1) for all $N$ states of $\bm{x}$, the probability density functions (PDFs), $P(\bm{x};t)$, can be obtained numerically. Since the focus of this research is the PDFs of Langevin equations with external forces, the Fokker-Planck equation (FPE), which gives the time evolution of the PDF, is used and is given by Eq. (2) Risken (1996),

[TABLE]

where $P(\bm{x},t)$ is the PDF of $\bm{x}$ at time $t$, and $\widehat{L}(\bm{x};t)$ is an FPE operator written as

[TABLE]

The time-dependent PDF can be obtained by solving Eq. (2).

II.1 Variational Superposed Gaussian Approximation (VSGA)

The VSGA is an adaptation of the multiple Gaussian wave packets approximation in quantum mechanics to Langevin equations Skodje and Truhlar (1984); Worth et al. (2004); Zoppe et al. (2005). It uses superposed multiple Gaussian distributions to approximate PDFs and obtain time-evolution equations for parameters of each Gaussian distribution with the variational principle Hasegawa (2015).

We first review procedures of the VSGA (for details of the method, please see Hasegawa (2015)). Starting from the FPE, Eq. (2), let $\Theta(\bm{x};t)$ be the time derivative of $P(\bm{x};t)$. Since $\Theta(\bm{x};t)dt+P(\bm{x};t)=P(\bm{x};t+dt)$, given $P(\bm{x};t)$ at time $t$, $P(\bm{x};t+dt)$ is specified when $\Theta(\bm{x};t)$ is determined. Thus we seek to obtain the optimal $\Theta(\bm{x};t)$ such that $P(\bm{x};t+dt)$ optimally satisfies the FPE (Eq. (2)) with given $P(\bm{x};t)$. The optimal $\Theta(\bm{x};t)$ should minimize

[TABLE]

When $P(\bm{x};t)$ satisfies Eq. (2), $R[\Theta]$ always vanishes. Ensuring the normalization condition ($\int P(\bm{x};t)d\bm{x}=1$) of the PDF yields a constraint condition on $\Theta(\bm{x};t)$ as $\int_{-\infty}^{\infty}\Theta(\bm{x};t)d\bm{x}=0$. To minimize Eq. (4) with the normalization condition equation, we introduce the Lagrange multiplier $\lambda(t)$ as follows:

[TABLE]

With the variation of $\delta\Theta(\bm{x};t)$ in Eq. (5), Eq. (6) should hold for optimal $\Theta(\bm{x};t)$,

[TABLE]

The approximate $P(\bm{x};t)$ obtained with the superposition of multiple Gaussian distributions is given as follows:

[TABLE]

where $\mathbf{A}_{m}$ is a positive definite $N\times N$ symmetric matrix, $\bm{b}_{m}$ is an $N$-dimensional column vector, and $r_{m}$ is the weight of the $m$th basis function. $N_{B}$ is the number of basis functions. In Eq. (7), a parameter set $\bm{\theta}(t)=(\theta_{1}(t),\theta_{2}(t),...,\theta_{K}(t))$ is composed of $\mathbf{A}_{m}$, $\bm{b}_{m}$, and $\bm{r}_{m}$. For an $N$-dimensional system and $N_{B}$ basis functions, the total number of parameters is

[TABLE]

Thus, Eq. (6) can be written as follows:

[TABLE]

We can uniquely obtain parameters of multiple Gaussian distributions by solving the $K$ constraints (Eq. (9)) and ensuring a normalization condition.

II.2 Augmented Variational Superposed Gaussian Approximation (A-VSGA)

As reported in Hasegawa (2015), the VSGA provides very accurate approximations for one- and two-dimensional driven systems. Still, this approach is limited in some respects; the integral in Eq. (9) can be computed analytically if the expectation of the drift term allowing the Gaussian integration, which enables an efficient calculation of time-dependent solutions. Typically, this condition is met when the drift term is given by polynomial functions because the expectation with respect to Gaussian distribution can be calculated analytically for arbitrary degrees of polynomials. However, the potentials of most non-linear models are not polynomial De Groot and Mazur (2013), and this makes the VSGA computationally expensive. Consequently, in the A-VSGA, we mitigate that restriction so that the approximation can handle rational polynomial functions.

We explain our method by reviewing a one-dimensional case with an additive noise term for simplicity. Essentially, the same calculation can be conducted for multivariate cases and for systems with multiplicative noise terms given by rational functions. We show that Eq. (9) can be solved when the drift term is a rational polynomial function. Suppose that the drift function is given by the rational function

[TABLE]

where $g(x)$ and $h(x)$ are polynomials. Substituting Eq. (10) into Eq. (2), we obtain

[TABLE]

Multiplying Eq. (11) by $h^{2}(x)$, we find

[TABLE]

When the time evolution of $P(x;t)$ satisfies Eq. (2) with the drift term (Eq. (10)), the right-hand side of Eq. (12) should vanish. Therefore, the optimal $\Theta(x;t)$ should minimize the following performance index:

[TABLE]

In a similar procedure to that of the VSGA, the Lagrange multiplier $\lambda(t)$ for the normalization constraint is introduced, and the following $3N_{B}$ equations are obtained (as Eq. (8) implies, the total number of parameters is $K=3N_{B}$ for one-dimensional systems):

[TABLE]

where

[TABLE]

We show that the A-VSGA can handle rational polynomial systems. Furthermore, it can be used to approximate non-linear systems with irrational potential functions. Our idea is that the drift terms can be approximated by a rational function via methods such as Taylor series expansion and the Padé approximant. Although Taylor expansion is an expression of functions as an infinite sum of monomials, it may not converge over the entire domain. Even in such a situation, the Padé approximant still works and often yields better approximation than truncating the Taylor series Stewart (1960); Fair (1964); Lorentz et al. (1996); Cohen (1991). A given function $V(x)$, can be approximated by the Padé approximant of order $[m/n]$ as follows:

[TABLE]

For systems with irrational potential, we first convert the drift terms to rational functions by using Padé approximants of a suitable order. From Eq. (14) and the normalization condition, we have $3N_{B}+1$ implicit differential equations. By numerically solving these differential equations (Mathematica 10 is used for the implementation), we can obtain accurate time-dependent solutions.

III Results

We apply the A-VSGA to a stochastic biochemical reaction and a soft bistable model. We also conduct Monte Carlo (MC) simulations to confirm the reliability of the A-VSGA. For MC simulations, we employ the Euler method with time resolution $\Delta t=0.001$, and the PDFs were constructed by repeating the stochastic simulations $100000$ times.

III.1 Gene Regulatory System

We first apply the A-VSGA to a model of genetic regulation with a positive auto-regulatory feedback loop. A flow diagram of the model is found in FIG. 1. The transcription factor activator (TF-A) gene incorporates a transcription factor responsive element (TF-RE). When phosphorylated TF-A homodimers bind to this element, the transcription is increased. The fraction of dimers phosphorylated is dependent on the activity of kinases whose activity can be regulated by external signals Liu and Jia (2004) (for more details of the model, please see Liu and Jia (2004)). The change rate of concentration of the TF-A dimer is described by the following Hill equation:

[TABLE]

where $k_{f}$ is the maximal rate of TF-A that activates transcription. TF-A is degraded at a rate of $k_{d}$ and synthesized at a rate of $R_{\mathrm{bas}}$. $K_{d}$ is the dissociation concentration of the TF-A dimer from TF-REs.

To simulate the stochastic effects of the synthesis process, we assume that white Gaussian noise, $\xi(t)$, is added to the reaction rate as $R_{\mathrm{bas}}\rightarrow R_{\mathrm{bas}}+\xi(t)$. Thus, Eq. (16) becomes the following Langevin equation:

[TABLE]

where $D$ is the noise intensity, $I(t)$ is an external driving signal, and $\xi(t)$ is white Gaussian noise. We set the parameters to $k_{f}=6$, $k_{d}=1$, $K_{d}=10$, and $R_{\mathrm{bas}}=0.4$, as is done in Liu and Jia (2004), to make the system bistable. The corresponding FPE operator $\widehat{L}(x;t)$ is given by Eq. (18),

[TABLE]

Substituting Eq. (18) into Eq. (14), we have $3N_{B}+1$ constraints including a normalization condition.

We first consider the stationary case, i.e., $I(t)$ is constant, since stationary PDFs can be calculated analytically for it. The analytic solution is given by Eq. (19) Er (1998),

[TABLE]

where $Z(D)=\int_{-\infty}^{\infty}\exp{\left[-U(x)/D\right]}dx$ and $U(x)=-\int\left[k_{f}x^{2}/\left(x^{2}+K_{d}\right)-k_{d}x+R_{\mathrm{bas}}\right]dx=-k_{f}x-R_{\mathrm{bas}}x+k_{d}x^{2}/2+k_{f}\sqrt{K_{d}}\arctan{\left(x/\sqrt{K_{d}}\right)}$.

FIG. 2 shows the stationary solutions obtained by the A-VSGA (blue line) and the analytic solutions (red line) obtained by Eq. (19) for (A) $D=1.0,I(t)=0,N_{B}=3$; (B) $D=1.0,I(t)=0.1,N_{B}=3$; and (C) $D=0.5,I(t)=0,N_{B}=3$. Note that there is a non-zero probability for $x<0$ that does not exist in reality since $x$ is the concentration of the TF-A monomer. Still, the purpose of this numerical experiment is to show that the A-VSGA can approximate PDFs of systems having rational functions. Therefore, the existence of non-zero probability in the negative region is not relevant for our purpose. For all the cases we have tested, the A-VSGA provides close PDFs as compared to those obtained analytically. In FIG. 2, the dashed lines describe each of the Gaussian basis functions that comprise the PDFs of the A-VSGA. Two bases each approximates one of the two peaks, while another basis is located at the center of the two peaks, reducing the error of the PDFs. Still, we see a small deviation from the analytical calculations in the A-VSGA results, which is the limitation of approximating PDFs with $N_{B}=3$ basis functions. For all parameter settings we have put in this experiment, the maximum feasible value of $N_{B}$ is 3. When the noise intensity approaches 0.1, it is difficult to find valid initial values. The nonorthogonality of multiple Gaussian distributions causes numerical instability, which also occurs in the original VSGA.

Stochastic processes with external forces are non-stationary processes. For systems driven by periodic forces, time-dependent solutions are often represented as a Fourier series expansion Jung (1993). However, Fourier series expansion cannot be adopted to solving Langevin equations with chaotic external forces. In the dynamic case, we add a chaotic input signal to the system to show that the A-VSGA can be used for such systems. Consider a Rössler oscillator:

[TABLE]

where we use $c_{1}=0.2,c_{2}=0.2$, and $c_{3}=5.7$. The chaotic input signal is defined as follows

[TABLE]

where $\alpha$ is the input strength and $\omega$ is the reciprocal of the time scale. FIG. 3 shows trajectories of aperiodic input signals $I(t)$. The parameter settings are $\alpha=0.02$ and $\omega=0.25$ (FIG. 3(A)) and $\alpha=0.02$ and $\omega=0.5$ (FIG. 3(B)).

FIG. 4 shows the PDFs as functions of $t$ and $x$, which are calculated by the A-VSGA with the parameter settings $D=1.0$, $\alpha=0.02$, and $\omega=0.25$ (FIG. 4(A)); $D=0.5$, $\alpha=0.02$, and $\omega=0.5$ (FIG. 4(B)); and $D=0.5$, $\alpha=0.02$, and $\omega=0.25$ (FIG. 4(C)). For non-linear systems driven by time-dependent external forces, there are no analytic solutions for Langevin equations. To investigate the accuracy of the $P(x;t)$ obtained by the A-VSGA, we performed MC simulations and compared the results of our method with that of MC simulations at a specific time. FIG. 5 shows the PDFs calculated by the A-VSGA (solid line), MC simulations (circles), and basis functions (dashed line) at time $t=100$. The parameter settings of FIG. 5 are the same as those for the data shown in FIG. 4. When $\omega=0.25$ (FIG. 5(A) and (C)), the PDFs of the A-VSGA agree very well with the results of MC simulations. The time-dependent solutions can be well approximated by 3 basis functions. In FIG. 5(B), the performance of the A-VSGA is not ideal. There is deviation between the curves around the peaks. The base near the left peak ($y=0$) is higher than the PDF obtained by the MC simulations, and the superposition of the two bases, near the right peak ($y=4.5$), is slightly lower than the solution obtained by the MC simulations. The results of our method would likely be more accurate in this case if more basis functions could be used.

III.2 Soft Bistable System

By using the Padé approximant, several non-linear functions that cannot be approximated by the Taylor series expansion can be approximated by rational polynomial functions. Then, we can use the A-VSGA to obtain their time-dependent solutions. Here, we consider a soft bistable system Longtin et al. (1994). The Langevin equation is

[TABLE]

where $I(t)$ is an input signal, $D$ is the noise intensity, and $\xi(t)$ is white Gaussian noise. The drift term is a hyperbolic function $\tanh{x}$ and not a polynomial.

The Padé approximant can provide a drift function that is very close to $\tanh(x)$, whereas the Taylor series expansion of the original function diverges due to a pole at $|x|=\pi/2$. When we use the Padé approximant to approximate the original function, we should pay attention to the following points: (1) the denominator of the approximant function should not have any poles, since poles lead the the A-VSGA to fail, (2) the Padé approximant of the drift term should approach $\infty$ for $x\to\pm\infty$, which may not be satisfied for some order $[m/n]$ (as FIG.6 shows, only odd order approximant of $\tanh{x}$ satisfy condition(2)), and (3) a high order $[m/n]$ may cause unavailability of the initial conditions and result in an increase in the calculation time. FIG.6 shows several Padé approximant of $\tanh(x)$ with different order $[m/n]$. In this case, when the order is odd, the approximant function satisfies both conditions (1) and (2). When the approximant order is $[5/5]$ or higher, there are more than 20 of the highest order of moment functions, which causes extremely long computation and numerical instability in the A-VSGA. Here, for the given function, we chose a Padé approximant of order $[3/3]$,

[TABLE]

Thus, the corresponding FPE operator $\widehat{L}(x;t)$ is given by,

[TABLE]

Substituting Eq. (26) into Eq. (14), we again have $3N_{B}+1$ constraints, including a normalization condition.

Again, we first consider a stationary solution of the soft bistable model where the analytical PDFs can be obtained by the method described in Er (1998). FIG. 7 shows the results of the our approach (blue line) and the analytic solutions (red line) obtained by Eq. (19) for $D=0.1$, $I(t)=0$, and $N_{B}=5$ (FIG. 7(A)); $D=0.1$, $I(t)=0.1$, and $N_{B}=5$ (FIG. 7(B)); and $D=0.5$, $I(t)=0$, and $N_{B}=3$ (FIG. 7(C)). For the analytical solutions, the potential function $U(x)=-\int(\tanh x-x/2)dx={x}^{2}/4-\ln\left(\cosh\left(x\right)\right)$ is inserted in Eq. (19). For $D\textgreater 0.5$, the PDFs can be approximated by 5 basis functions (dashed line). Two bases each are located near the deterministic stable steady states ($y=\pm 2$), and one is near the deterministic unstable steady state ($y=0$). The maximum feasible value of $N_{B}$ is clearly larger than that in the previous experiment. The stationary PDFs obtained by the A-VSGA are accurate. A large number of basis functions significantly improves the accuracy of the results. The A-VSGA also shows very good agreement with the analytical solutions in panel FIG. 7(C), though only 3 basis functions are used. In this experiment, the highest order of moment function in Eq. (14) is the same as that in the first experiment, the gene regulatory system. However, here the results are more accurate. The only difference between the two experiments is that the complexity of Eq. (14) of the soft bistable system is lower than that of the gene regulatory system. So, we suggest that a complex Eq. (14) may increase the calculation difficulty, thus resulting in extensive temporal cost.

Next we study the time-dependent solutions of a driven system. We add a chaotic input signal (defined by Eq. (23)) to the soft bistable system. FIG. 8 shows the PDFs of the soft bistable system driven by the external input signal $I(t)$ with the following parameter settings: $D=0.2$, $\alpha=0.02$, and $\omega=0.5$ (FIG. 8(A)); $D=0.2$, $\alpha=0.02$, and $\omega=0.25$ (FIG. 8(B)); and $D=0.5$, $\alpha=0.02$, and $\omega=0.25$ (FIG. 8(C)). FIG. 9 shows the time-dependent solutions, $P(x;t)$, at $t=100$, which are calculated by the A-VSGA (solid line) and MC simulations (circles) with the same parameter settings used for the PDFs of the soft bistable system (FIG.8). When noise intensity $D=0.2$, the PDFs can be approximated by 4 Gaussian distributions (FIG. 9(A) and (B)). Near each of the two peaks ($y=\pm 2$) is a base, and at the center of each of the two peaks ($y=0$) there is a base. For all cases we have tested, the results of the A-VSGA are in excellent agreement with results of the MC simulations. In FIG. 9 (C), the PDF is a superposition of 3 Gaussian distributions. This verifies the reliability of the A-VSGA with respect to the PDFs at a specified time in dynamic case.

IV Conclusion

We propose the A-VSGA for the time-dependent solutions of Langevin equations with rational or several other non-linear drift terms. Our approach is applied to two systems to show its reliability. Compared with the MC simulations, the A-VSGA is an efficient method to calculate time-dependent solutions without relying on stochastic approaches. The results of the soft bistable system in Sec. III are satisfactory and agree well with the PDFs from the MC simulations.

Our approach uses multiple Gaussian distributions to approximate PDFs, which is the same as is done in the VSGA. A large number of basis functions that are not orthogonal prevents the calculation of the time evolution of the parameters. It can be seen from Sec. II that soving Eq. (14) in the A-VSGA is more complicated than solving Eq. (9) in the VSGA. The complexity of Eq. (14) increases the computation difficulty. Considering the two experiments carried out in Sec. III, the maximum feasible value of $N_{B}$ of the soft bistable system is five, whereas no more than three basis functions are used in the gene regulatory system. The reason is that Eq. (14) in the gene regulatory system is more complex than that of the soft bistable system. Moreover, an overly complex calculation in the gene regulatory system results in no valid initial values for the A-VSGA when the noise intensity is relatively small. Thus, the complexity of Eq. (14) is related to the time cost and difficulty of the calculation.

The A-VSGA can be used to obtain time-dependent solutions of Langevin equations whose drift terms are expressed by rational function or other forms such as a hyperbolic function. We have shown that the A-VSGA can effectively provide accurate approximation in non-equilibrium systems. This approach can be used for analyzing certain real-world problems such as enzymatic reactions, cell cycles, gene regulation, and a variety of ligand-receptor interaction problems.