Construction of quasi-potentials for stochastic dynamical systems: an optimization approach
CConstruction of quasi-potentials for stochastic dynamical systems: an optimizationapproach
R. D. Brackston, ∗ A. Wynn, and M. P. H. Stumpf Department of Life Sciences, Imperial College London, SW7 2AZ, United Kingdom Department of Aeronautics, Imperial College London, SW7 2AZ, United Kingdom (Dated: May 21, 2018)The construction of effective and informative landscapes for stochastic dynamical systems hasproven a long-standing and complex problem. In many situations, the dynamics may be described bya Langevin equation while constructing a landscape comes down to obtaining the quasi-potential, ascalar function that quantifies the likelihood of reaching each point in the state-space. In this work weprovide a novel method for constructing such landscapes by extending a tool from control theory: theSum-of-Squares method for generating Lyapunov functions. Applicable to any system described bypolynomials, this method provides an analytical polynomial expression for the potential landscape,in which the coefficients of the polynomial are obtained via a convex optimization problem. Theresulting landscapes are based upon a decomposition of the deterministic dynamics of the originalsystem, formed in terms of the gradient of the potential and a remaining “curl” component. Bysatisfying the condition that the inner product of the gradient of the potential and the remainingdynamics is everywhere negative, our derived landscapes provide both upper and lower boundson the true quasi-potential; these bounds becoming tight if the decomposition is orthogonal. Themethod is demonstrated to correctly compute the quasi-potential for high-dimensional linear systemsand also for a number of nonlinear examples.
I. INTRODUCTION
Multi-dimensional nonlinear dynamical systems canexhibit complex behavior including limit cycles, strangeattractors and multiple fixed points. When also drivenby stochastic perturbations, such systems often explorea range of conditions and may exhibit behavior such asrandom switching between attractors and stochastic reso-nance [1]. The mathematics of such systems has its rootsin the description of Brownian motion, which motivatedthe study of stochastic differential equations, generallywritten in the form,d x = f ( x ) d t + g ( x ) d W t . (1)Here the deterministic component of the dynamics (thedrift) is described by f ( x ), while the stochastic com-ponent (diffusion) is given by g ( x ) d W t , where d W t de-scribes the increment of a Wiener process [2].A popular description of such systems is in terms of anenergy landscape, in which the dynamics are describedas a ball moving in a potential basin. In some contexts,the landscape may simply give an intuitive description ofthe dynamics [3, 4], while in others it may represent atrue energy function [5]. In the context of developmen-tal biology, Waddington’s epigenetic landscapes providea popular analogy for stem cell development [6, 7], eventhough it has proven impossible to define a true energyfunction for general developmental processes.If correctly formed, landscapes may offer a quantita-tive analysis and interrogation tool beyond a merely phe-nomenological description, even in cases where the free ∗ [email protected] energy of the system may not be defined. In particular,a so-called quasi-potential landscape may be defined interms of the Freidlin-Wentzell action functional [8]. Thislandscape provides quantitatively a measure of the rel-ative stability of different states and the most probablepaths between them. The concept of landscapes definedby the action has received particular recent attention[9], as has the development of methods to calculate theso-called Minimum Action Path (MAP) [10–12]. Otherrecent methods have sought to evaluate quasi-potentiallandscapes with a focus on the action functional. For lin-ear systems an analytic expression for the quasi-potentialwas obtained by [13], while numerical methods applicableto nonlinear systems have been given in [14–16]. Whilesuch methods may offer a numerical solution over a dis-cretized space, in this work we present a method thatgenerates an analytical solution, and furthermore is ap-plicable to both linear and nonlinear systems.Given the motivation behind landscape descriptions, itis unsurprising that considerable effort has been appliedto find other methods for landscape evaluation basedeither upon experimental data or mathematical mod-els [17, 18]. In biological applications, the most popu-lar and readily applied method is that based upon thesteady-state probability distribution [19]. This approachis based upon the Fokker-Planck equation, which justi-fies computing a potential as U ( x ) ∝ − ln( P S ( x )). Inpractice the steady-state distribution is either obtainedby solving the Fokker-Planck equation directly [20–22],or found from extensive simulations of the correspondingSDE [23–25]. While this method is easy to implement, itmay be impractical for higher-dimensional systems andthere is generally no guarantee that the derived landscaperelates directly to the fixed points of the deterministicsystem. a r X i v : . [ m a t h . O C ] M a y Irrespective of how a landscape is obtained, its exis-tence implies a decomposition of the deterministic vectorfield into two components, one of which is given by thegradient of the landscape while the other accounts for theremainder of the dynamics in f ( x ): f ( x ) = −∇ U ( x ) + f U ( x ) . (2)The construction of the landscape may therefore also beperformed by considering the properties of this decom-position directly. As discussed by [8] and advocated by[26], one particularly interesting case is that in whichthe gradient and remainder are everywhere orthogonal,implying an independence between gradient-based androtational dynamics. While this property may be sat-isfied in some cases, it is important to note that it re-mains unclear if this is always possible. Motivated bythis, we develop here a method to form landscapes satis-fying the slightly more general sub-orthogonal conditionthat ∇ U ( x ) · f U ( x ) ≤
0, of which the orthogonal decom-position is a special case.As we shall discuss below, the sub-orthogonal decom-position has the desirable property that under certainconditions it is directly related to the quasi-potential.One of the other desirable properties of the landscape,and one that is satisfied by the sub-orthogonal decompo-sition, is that it should be a Lyapunov function [27] forthe deterministic dynamics. In this context, Lyapunovfunctions are often used to prove the asymptotic stabil-ity of the fixed points f ( x ) = 0. However, they alsohave the useful property that they provide the basinsof attraction for the deterministic dynamics. There istherefore a strong equivalence between quasi-potentiallandscapes and Lyapunov functions, suggesting that pre-existing methods for generating such functions may alsobe used to obtain the quasi-potential. This is the ap-proach that we take here.In this paper, we develop a method that performs asub-orthogonal decomposition of the deterministic dy-namics f ( x ), by utilizing a popular method for the con-struction of Lyapunov functions: the Sum-of-Squares(SOS) method [28]. We will firstly provide some back-ground on the properties of the sub-orthogonal decompo-sition in section II before providing details of this compu-tational method in section III. We will then examine itsapplication to a series of examples. The first of these ex-amples will be linear systems in section IV for which ref-erence analytical solutions exist, followed by applicationto nonlinear systems with various properties in section V.Conclusions will finally be given in section VI. II. PROPERTIES OF THE SUB-ORTHOGONALDECOMPOSITION
In this work we consider the case in which the deter-ministic part of (1) is decomposed as in (2), and thisdecomposition has the property that, ∇ U ( x ) · f U ( x ) ≤ . (3) When this relation becomes an equality, the decompo-sition can be said to be orthogonal, while otherwise wewill refer to the decomposition as sub-orthogonal . As weshall now discuss, such a decomposition of the vector fieldhas a useful relation to the quasi-potential and further-more provides a Lyapunov function for the deterministicdynamics. A. The quasi-potential
Let us now restrict our analysis to systems in whichthe diffusion tensor g ( x ) is a uniform diagonal matrix,describing purely additive noise,d x = f ( x ) d t + σI n d W t . (4)Here, σ is a small constant and I n is the n × n identitymatrix. Given such a system, the action associated witha path x ( t ) = ϕ is given by the following integral, S ( ϕ ) = 1 σ (cid:90) T | ˙ ϕ − f ( ϕ ) | d t. (5)Equation (5) is useful because in the limit, σ → P ( ϕ ) ∝ e − S ( ϕ ) . (6)Given the properties of the action, one can define aquasi-potential Q a with respect to a fixed point a as, Q a ( x ) = inf ϕ,T [ S ( ϕ ) | ϕ (0) = a, ϕ ( T ) = x ] . (7)This quasi-potential therefore gives the action associatedwith the minimum action path to each point within thebasin of attraction of the fixed point, and in turn maybe related to the probability distribution over the statespace.Given the sub-orthogonality condition of (3), thequasi-potential may be evaluated as, Q a ( x ) = 1 σ inf ϕ,T (cid:34)(cid:90) T | ˙ x + ∇ U − f U | d t (cid:35) = 1 σ inf ϕ,T (cid:34)(cid:90) T | ˙ x − ∇ U − f U | +4 ˙ x · ∇ U − f U · ∇ U d t (cid:35) = 4 σ ( U ( x ) − U ( a ))+1 σ inf ϕ,T (cid:90) T | ˙ x − ∇ U − f U | (cid:124) (cid:123)(cid:122) (cid:125) (i) − f U · ∇ U (cid:124) (cid:123)(cid:122) (cid:125) (ii) d t In the orthogonal case, term (ii) is equal to zero, whilethe infimum of term (i) will approach zero for a patharbitrarily close to that satisfying the ODE ˙ x = ∇ U + f U [29]. For the sub-orthogonal case the net contributionsof both terms will be positive. The potential U thereforeprovides a lower bound for the quasi-potential accordingto, Q a ( x ) ≥ σ ( U ( x ) − U ( a )) . (8)In the case of the orthogonal decomposition this inequal-ity becomes an equality. B. Lyapunov functions
Lyapunov functions serve an important purpose in thefield of non-linear dynamical systems and feedback con-trol, in proving asymptotic stability of a fixed point.While the proof of asymptotic stability requires specificconditions to be satisfied, here we use a slightly looserdefinition that allows the Lyapunov function to havewider interpretation. For a system defined by a set ofODEs ˙ x = f ( x ) a Lyapunov function U ( x ) : R n → R , isone that satisfies the following constraints: U ( x ) ≥ , (9a) ∇ U ( x ) · f ( x ) ≤ . (9b)These constraints impose the requirements that U is ev-erywhere positive and that U decreases along trajecto-ries of f ( x ). The key implication of these constraints isthat the Lyapunov function correctly captures the basinsof attraction for all the stable fixed points: a region ofthe state space A around a fixed point a such that if x (0) ∈ A , x ( t ) → a as t → ∞ . Such a function thereforeprovides an accurate qualitative description of a land-scape underlying the dynamics.Given the sub-orthogonal decomposition property of(3), we may rewrite constraint (9b) as, ∇ U · f = ∇ U · ( −∇ U + f U )= −|∇ U | + ∇ U · f U ≤ . (10)Sub-orthogonality of the decomposition is therefore a suf-ficient condition for the potential U to be a Lyapunovfunction for the deterministic dynamics described by f .We therefore use this condition in our optimization ap-proach discussed below. III. THE OPTIMIZATION METHOD
In this work we develop a method for landscape for-mation based on the sum-of-squares method for formingLyapunov functions. We will therefore first give a briefoverview of the standard algorithm through which SOSis used to find Lyapunov functions, before detailing the key changes required in our algorithm to achieve a sub-orthogonal decomposition. The algorithm itself has beenimplemented in the programming language Julia, builtupon the JuMP package for mathematical optimization[30] and in
Matlab using the SOStools package [31].
A. The standard SOS method
Finding a Lyapunov function for nonlinear dynamicalsystems is generally a difficult problem[27, 32], and isknown to be NP-hard in the case where U is a polynomial[33]: i.e. given a possible function U , even checking thatit satisfies the requirements may be very computationallyexpensive. The basis of the SOS method is that polyno-mial functions formed as the sum of the squares of lowerorder polynomials are guaranteed to be positive-definite,and furthermore can be found via efficient optimizationmethods. While the requirement that U ( x ) is SOS isstricter than that of positive-definiteness, it makes theproblem computationally tractable.The standard use of SOS in forming a Lyapunov func-tion involves solving the following:find a feasible [ c j ] mj =0 subject to U = c + m (cid:88) j =1 c j p j (11a) U ≥ (cid:15) n (cid:88) i =1 x i , (11b) ∇ U · f ≤ . (11c)Here the p j are a suitable set of m monomial terms (e.g. x x ) and the c j are real coefficients. The parameter (cid:15) is a positive constant, typically of order one. The tasktherefore involves a pure feasibility problem, and comesdown to finding any U composed of the monomials spec-ified in (11a), subject to the inequality constraints (11b),(11c) [34]. There are therefore infinitely many possi-ble solutions, if f ( x ) describes a stable system. As wewill discuss below, we modify each of the steps of prob-lem (11) to produce a convex optimization problem thatyields a unique result. B. Towards orthogonality
While the SOS method is able to generate Lyapunovfunctions for general n -dimensional dynamical systems,such functions are not unique and will generally notsatisfy the sub-orthogonality requirement of (3). Asnoted by several previous authors [14, 26], an orthogo-nality requirement results in a Hamilton-Jacobi equation ∇ U · f + ||∇ U || = 0 which is a nonlinear constraint on U ( x ). Such constraints cannot be directly implementedinto the SOS program. In order to implement such aconstraint in a linear manner, consider the matrix M U := (cid:20) −∇ U · f ∇ U (cid:62) ∇ U I n (cid:21) ∈ R ( n +1) × ( n +1) . (12)A Schur Complement argument [35] implies that M U (cid:23) x ∈ R n if and only if − ∇ U · f − ∇ U (cid:62) I n ∇ U ≥ , (13)which is equivalent to ∇ U · f U ≤
0, where f U is as definedin (2). This in turn implies that ∇ U · f ≤
0, whichsatisfies the key requirement for a Lyapunov function andconstraint (11c). With the matrix inequality constraintin place, there will still be an infinite number of possible U , most of which will not come close to orthogonality.In order to find the unique U that minimizes |∇ U · f U | ,we find it necessary to make U as steep as possible, inthis way maximizing the contribution of −∇ U to f . Thisis achieved by maximizing an appropriately chosen lowerbound, B ( x ) = (cid:88) i (cid:15) i b i ( x ) , (14)where each b i is a monomial in x and the (cid:15) i are realcoefficients. The choice of these monomials is importantand the method for choosing them is described belowin section III C. The complete optimization procedure istherefore as follows:maximize [ c j ] , [ (cid:15) i ] (cid:88) i (cid:15) i (15a)subject to U = c + m (cid:88) j =1 c j p j (15b) U ≥ (cid:88) i (cid:15) i b i , (15c) M U (cid:23) . (15d)For the case that f ( x ) is linear, this method can beproven to yield the correct result, as detailed in ap-pendix A. C. Choosing the monomial basis and lower bound
The first key choice that must be made is that of themonomial basis from which U is constructed, i.e. the p j in (15b). In the typical application of SOS to form Lya-punov functions, a candidate function is formed from afull collection of the monomials in x up to a given (even)degree d [36]. For a given n -dimensional system, this re-sults in m = (cid:0) n + dd (cid:1) monomials forming the polynomialfunction U [37]. For example, for a 4th order polynomialand 4-dimensional state space this results in 40 mono-mials. Since the computational cost of the optimizationincreases with m , it is advantageous to reduce the ba-sis if possible, by exploiting particular properties of thefunction we are trying to obtain.
1. A minimal basis
We can initially provide an upper bound on the degree d using the sub-orthogonality condition (3), which maybe rewritten as, ∇ U · f + |∇ U | ≤ . (16)Since |∇ U | ≥
0, the above inequality can only hold ifthe degree of ∇ V · f is at least as big as the degree of |∇ U | . If d is the degree of U and e the degree of f , then d + e − ≥ d −
2, meaning that, d ≤ e + 1 . (17)While (17) provides an upper bound on the degree of U , we may further refine the basis (15b), motivated bythe properties in the pure gradient case, in which f i = − ∂U/∂x i . The minimal basis for such a U is thereforeprovided by [ p j ] mj =1 = M (cid:34)(cid:88) i f i x i (cid:35) , (18)where the operator M extracts a vector of monomialsfrom a polynomial.As an example, consider the vector field defined by, f ( x ) = x − x − x − x . (19)Here the highest order in f is three so (17) implies that wewill likely require a d = 4 order polynomial. Given thatthe system has n = 3 dimensions, the standard methodleads to 35 monomials. By contrast, a minimal basiswould only include five terms, [ x , x , x , x ,
2. The lower bound
Given a monomial basis for U , it is next necessary tochoose the monomials b i in the lower bound (14). Thischoice is best illustrated via the following simple systemdescribed by, f ( x ) = x − x + β = − dd x (cid:20) x − x − βx + C (cid:21) . (20)Given that (20) is one-dimensional and therefore maybe written as a pure gradient system, the minimal ba-sis argument of (18) correctly ascertains that U ( x ) willbe composed only of monomials [ x , x , x, b ( x ) = x q , where q isan even integer. The inequality constraint (15c) thereforebecomes: c x + c x + c x + c ≥ (cid:15)x q , (21)where the choice of q is of critical importance. In partic-ular it must satisfy two particular properties: FIG. 1. A schematic of the lower bounding method. (a) Thepotential function in the region of the fixed points and (b) acomparison between the potential and lower bound. a. Sufficient pressure
Suppose that we take q = 2,now for any c , c > c may remain arbitrarily small and the problem willnot have converged. The lower bound is therefore unableto exert sufficient “upward pressure” on U . b. Feasibility Alternatively consider q = 6. In thiscase it is impossible to find coefficients such that U > (cid:15)x ,since as x → ∞ , c x < (cid:15)x for any c , (cid:15) >
0. Theprogram will therefore be infeasible.The correct choice for q in this case is q = 4, i.e. thehighest degree in the basis for U . An example of thislower bound in practice is displayed in figure 1. Thesame principles of sufficient upward pressure and feasi-bility also apply to higher dimensional systems, such thatthe lower bound is chosen to consist of the highest degreesingle and mixed monomials in the basis for U .
3. Extending the basis
A correctly chosen “minimal” basis will naturally workfor the case that f ( x ) is a pure gradient system, and may sometimes work in other cases [38]. However more gener-ally, our constraints may require that ∇ U contains mono-mials beyond (18), as will be illustrated in later examples.Nonetheless, in addition to considerations of computa-tional cost, the sensitivity of the method to the choiceof lower bound means that we cannot simply choose allmonomials up to degree d , as is done in the standard SOSmethod.Consider again the example given in (20). Based onthe above discussion it is clear that we could not have in-cluded an additional higher-order monomial in the basisfor U . For example had we included a x term, we wouldlogically have also chosen q = 6 for the lower bound.Since the actual potential need not include a x term,this would make the program infeasible. The principlesof sufficient pressure and feasibility in the choice of thelower bound, therefore guide which monomial terms arepermitted in the basis for U . The procedure for deter-mining the basis is therefore as follows:1. The highest order monomial in any individual x i should be that obtained from the minimal basis.2. Other mixed monomials are required, providedthat:(a) the total order is not more than the highesttotal order in the minimal basis,(b) the individual order in each x i is not morethan the highest individual order monomialin the minimal basis. D. Iterative improvement
The process outlined above performs well for many sys-tems, but may require further improvement to push theobtained landscape closer to orthogonality. Given an ini-tial Lyapunov function U such as that found by imple-menting (15), one can attempt to iteratively improve thesolution, generating at each step a second Lyapunov func-tion U composed of the same monomial basis. Each it-eration involves solving the following optimization prob-lem:minimize [ c j ] ,α,(cid:15) α (22a)subject to U = c + m (cid:88) j =1 c j p j (22b) U ≥ (cid:15) (cid:88) i x i , (22c) M U (cid:23) , (22d) ∇ U · ( f + 2 ∇ U ) ≥ α ( f · ∇ U ) + (1 + α ) ||∇ U || , (22e) α, (cid:15) > . (22f)The key to this optimization is that by using the ini-tial guess, U , the optimization is constructed to be lin-ear in the new improved Lyapunov function, U . Con-straints (22c) and (22d), simply mirror those in optimiza-tion (15), namely positive definiteness of the Lyapunovfunction and the sub-orthogonality condition. Con-straint (22e) is chosen such that for α = 1 equality isguaranteed if U = U , while for α < ∇ U · f U ≤∇ U · f U . This result is proven in Appendix B, howeverthe key point is that the optimization has a guaranteedfeasible solution with α = 1 and U = U , while if asmaller α is obtained then the new Lyapunov function iscloser to orthogonality.The complete method for obtaining the potential land-scape therefore consists of first applying procedure (15),followed by a number of iterations of procedure (22).We generally find that only one or two iterations of thismethod are required for satisfactory convergence. Thenumber of iterations may therefore either be pre-specifiedor determined based upon a convergence criterion. IV. THE LINEAR CASE
A special case of the orthogonal decomposition may beconsidered when the deterministic system is linear [14].While this may seem a very limiting test case, the linearsituation can exemplify some of the key scenarios andlimitations that also apply in the nonlinear case, whileremaining tractable to analysis.
A. Properties of linear systems
Linear dynamics may be written in terms of a matrixmultiplication as, ˙ x = f ( x ) = A x , (23)where the matrix A ∈ R n × n . The construction of thequasi-potential then comes down to decomposing the ma-trix A into two parts as A = A g + A c , (24)where the gradient matrix A g is a symmetric matrix thatgives the potential according to U ( x ) = − x (cid:62) A g x . (25)Because A g is symmetric, it has purely real eigenvalues,while if the decomposition is orthogonal, the remaindermatrix A c has purely imaginary eigenvalues and is suchthat the product A g A c is an antisymmetric matrix.In such linear cases it has been demonstrated that anorthogonal decomposition does always exist [39], and fur-thermore may be provided by an analytical expression[13] as, A g = 12 (cid:18)(cid:90) ∞ e At e A ∗ t d t (cid:19) − . (26) If the matrix A is normal ( AA ∗ = A ∗ A ), then this ex-pression simplifies to, A g = 12 ( A + A ∗ ) x . (27)Since these expressions provide solutions against whichour method can be tested, linear systems give a perfecttest case for our optimization-based approach which doesnot rely on any prior knowledge of the properties of thesystem. B. A three-dimensional example
We consider now the test-case of a three-dimensionallinear system defined by, f ( x ) = A x = − . . . . − . . . − . − . x . (28)Implementation of methods (15) and (22) provides thefollowing decomposition for A : A = − .
01 0 .
14 0 . . − . − . . − . − . (cid:124) (cid:123)(cid:122) (cid:125) A g + . − .
14 0 . − .
14 0 .
05 3 . . − . − . (cid:124) (cid:123)(cid:122) (cid:125) A c , (29)giving the quasi-potential as, U ( x ) = 2 . x + 0 . x + 0 . x − . x x − . x x + 0 . x x . (30)As can readily be observed, the matrix A g is symmetricand can easily be verified to have purely real (negative)eigenvalues. The matrix A c is not itself either symmet-ric or antisymmetric but can be verified to have purelyimaginary eigenvalues. The dynamics associated with A c are therefore those of pure oscillation, and evolve on levelsets of the potential U .It is noteworthy that for the system defined by A , thereis no direct connection between x and x . One maytherefore naively expect such a direct connection to alsobe absent in the potential. Nevertheless, in order to pro-vide an orthogonal decomposition, the potential matrix A g does include such terms that are then negated in A c .This confirms that a truly minimal basis is indeed in-sufficient for forming the quasi-potential, as discussed insection III C 3. FIG. 2. Scaling of the computational cost with system size.Tests were performed on a basic desktop PC with 16 Gb ramand an Intel Core i3 processor.
C. Scaling with system size
While the method may be verified to work for thisthree-dimensional example, it is also useful to assess theaccuracy and efficiency of the algorithm for higher di-mensional systems. We do this by randomly generatingreal, negative-definite matrices of varying size n = 2 : 10.In all such cases the method is able to perform a de-composition that satisfies the orthogonality constraint.Results for the computational cost are shown in figure 2,displayed on a logarithmic scale. It is evident that thecost increases considerably with system size, most likelyin a factorial fashion. For small systems the algorithmconverges in less than a second, while for larger systemsit may take up to an hour.It is worth noting here that some of the additionalcomputational cost associated with larger systems arisesfrom a greater number of times that the iteration pro-cedure (22) must be applied, given a desired level of or-thogonality. For low-dimensional problems we find thatsometimes no iterations are required, while for systems ofdimension n = 10, up to five iterations may be necessaryto achieve the same quality of result. V. THE NONLINEAR CASE
Having studied the ability of our method to obtain thedecomposition in linear cases, we now move on to fullynon-linear examples for which analytical solutions do notalways exist.
A. A nonlinear multistable system
Our first example is the quartic system from [26] withfour attractors and a known orthogonal decomposition.
TABLE I. Coefficients obtained for the Maier-Stein model forthree different parameter combinations. µ γ c c c c The dynamics are defined by, f ( x ) = (cid:20) − x − x + 9 x − x − x + 2 x + 11 x − x (cid:21) , (31)for which the true potential is given by, U ( x ) = 0 . x + x ) − x + x ) + x x + x . (32)For such a problem our method finds the potential towithin 5 significant figures within a few seconds. Becausethe output of the algorithm is a symbolic expression forthe potential, the actual value of U may then be eval-uated at any required points in the space with minimalfurther effort. This is in contrast of course, to methodsthat solve the associated PDE over a discretized grid,since evaluation outside of the solved area requires sig-nificant further computation. B. The Maier-Stein model
We next apply the method to the widely studied Maier-Stein model of [40], defined by, f ( x ) = (cid:20) x − x − γx x − µ (1 + x ) x (cid:21) . (33)This model has the property that when γ = µ , the dy-namics can be expressed in terms of a pure potential [41] f ( x ) = −∇ U , where, U ( x ) = +0 . (cid:124) (cid:123)(cid:122) (cid:125) c x +0 . µ (cid:124) (cid:123)(cid:122) (cid:125) c x x − . (cid:124) (cid:123)(cid:122) (cid:125) c x +0 . µ (cid:124) (cid:123)(cid:122) (cid:125) c x . (34)For cases in which γ (cid:54) = µ , the system also has a non-gradient component.We apply our method for three different parametercombinations, obtaining in each case a polynomial ex-pression for the potential with the same non-zero terms.The coefficients for these terms are displayed in table I.For the non-gradient cases, the decompositions are, asdesired, sub-orthogonal. C. Bounds on the quasi-potential
A useful property of the (sub-)orthogonal decomposi-tion is that it may be used to obtain predicted MAPs. Asdiscussed in section II, MAPs are those paths that min-imize the action functional and provide a definition ofthe true quasi-potential. In general these paths must befound via a costly optimization that must be performedfor every point in the state-space. For a system that sat-isfies the orthogonal decomposition however, the MAPfrom a fixed point x o to another point x e within thebasin of attraction, is one that follows,˙ x = ∇ U + f U . (35)Such a path may therefore readily be obtained by a sim-ulation starting at x e and following ˙ x = −∇ U − f U .While we expect paths described by (35) to match ex-actly in the orthogonal case, for sub-orthogonal cases wehope that there may still be close agreement, depend-ing on the degree of orthogonality. We therefore com-pare the predicted and exact paths, evaluating the trueMAPs via our own implementation of the optimizationapproach detailed in [12]. A comparison of paths fromthe fixed point x o = ( − . , .
0) is given in figure 3. Inthe gradient case, the paths match exactly, as expected.For the first non-gradient cases, there is clearly some dis-agreement, although the sense of curvature of the pathsis generally the same for the predicted and true MAPs.While the agreement of the paths may seem quanti-tatively poor, these paths may still be used to providean estimate of the true quasi-potential. Given that thequasi-potential is defined in terms of the infimum of theaction over all possible paths, the action of any path thatwe choose can provide an upper bound to the true quasi-potential, evaluated via the geometric minimum actiomethod [10]. While most choices of possible path will givean action much greater than that of the MAP, the SOSpredicted path may be expected to give a tighter bound,owing to the qualitative similarity with the true MAP. Inaddition to this upper bound on the true quasi-potential,the potential from the decomposition itself provides alower bound, as discussed in section II.A quantitative comparison between the true quasi-potential Q , the potential from the decomposition U and that from the SOS predicted path S is given in fig-ure 4. In the gradient case there is almost exact agree-ment between all three methods, demonstrating that theobtained decomposition provides both the true quasi-potential and predicts the correct paths. In the firstnon-gradient case, the upper bound remains remarkablytight, despite the large differences between the true andpredicted paths. For the second non-gradient case thebounds are even tighter, with an exact match along the x -axis. In all cases, each of U and S can be seen to indeedbe lower and upper bounds respectively. VI. CONCLUSIONS
Evaluating the quasi-potential for linear and nonlinearSDEs is a challenging problem, yet one with significantinterest and motivation. In this work we have provided a novel method for the calculation of the quasi-potentialbased upon the Sum of Squares method for constructingLyapunov functions. Our method is applicable to sys-tems for which the governing equations are polynomialand involves solving an optimization over the coefficientsof a polynomial potential function.The construction of an informative landscape is moti-vated by three key requirements. Firstly, we would likethe potential to correctly capture the basins of attrac-tion for the deterministic system. Such a requirement isequivalent to the potential being a Lyapunov function,and is therefore naturally achieved by our method. Sec-ondly, it is desirable to have an estimate of the most prob-able transition trajectories between basins of attraction;the so-called minimum action paths. For cases permit-ting an orthogonal decomposition of the dynamics, thepaths may readily be obtained from the two vector-fieldcomponents. For cases in which the decomposition issub-orthogonal, the obtained paths may provide a more-or-less accurate approximation, suitable for use as an ini-tial guess in an optimization routine. Finally, we maywish for the potential function to be a quasi-potential forthe system, accurately describing the transition probabil-ities for situations of vanishingly small noise. Again, forsystems permitting an orthogonal decomposition of thedynamics our method calculates exactly such a quasi-potential, while for sub-orthogonal cases the potentialmay be used to provide both a lower and upper bound.The first key limitation of the method is in its applica-bility to only polynomial systems. However, such systemsare commonplace in e.g. mass action models of chemi-cal kinetics and linear models for dynamical systems. Aclosely related limitation is in the ability of the method toexpress the potential itself in terms of a polynomial. It isplausible that in some cases even if the governing equa-tions are polynomial, a potential satisfying the normaldecomposition must be expressed in some expanded ba-sis beyond monomial terms. Regardless, our polynomialsub-orthogonal potential still provides useful bounds.A second key limitation is that the obtained quasi-potential is only valid for systems with additive noise, inwhich the noise tensor is equal to the identity matrix.Yet this is a common approximation, especially whenthe magnitude of the noise tends to zero, and further-more can always be achieved for linear systems via acoordinate transform. If multiplicative noise is presentbut may expressed as a polynomial, it is possible that insome cases this could be incorporated into the algorithm,however this is beyond the scope of this study.While we have provided a method that generates ananalytical expression for the quasi-potential for a partic-ular subset of SDE systems, it remains unclear if such afeat can be achieved in more general cases. Ultimately,we hope that the method provides another useful tool forthe interpretation and analysis of stochastic dynamicalsystems, and may pave the way for more general meth-ods to obtain the quasi-potential in a symbolic manner.
FIG. 3. A comparison of the SOS estimated and true MAPs. Paths start at the fixed point x o = ( − . , .
0) and end at thepoints x e = ( − . , .
0) and x e = ( − . , . U is the potential formed from the sub-orthogonal decomposition, Q is the true quasi-potential and S is the action of the predicted minimum action path. ACKNOWLEDGMENTS
We are grateful to Dr David Schnoerr for critical read-ing of the manuscript. This work was supported by BB-SRC grant BB/N003608/1.
CODE AVAILABILITY
The algorithms discussed in this work have been im-plemented in both Julia and
Matlab , freely acces-sible at https://github.com/rdbrackston/SDEtools and https://github.com/rdbrackston/normalSoS re- spectively.
Appendix A: Optimality in the linear case
In section III B, we presented a method by which wemaximize a lower bound for the quasi-potential U , at-tempting to make the potential as steep as possible. Herewe justify this method by proving that in the linear case,the matrix defining a normal decomposition is also thatwhich provides the largest potential, defined in a suitableway.0 Lemma 1
Consider a full rank real matrix A ∈ R n × n whose eigenvalues have negative real part, and U ( x ) = x (cid:62) P x for some P = P (cid:62) ∈ R n × n .Then the decomposition, A x = −∇ U ( x ) + f U ( x ) , x ∈ R n , holds with f U ( x ) = ( A + P ) x .Suppose further that P solves the optimization problem:maximize P tr ( P ) (A1a) subject to P = P (cid:62) , (A1b) P ( A + P ) (cid:22) . (A1c) Then
P A + A (cid:62) P + 2 P = 0 . Hence, ∇ U · f U = 0 .Proof . First note that (A1) is feasible since A is stable.Furthermore, for any P satisfying P ( A + P ) ≤
0, we have, C := 2 P + P A + A (cid:62) P ≤ . Hence, the Ricatti equation,2 P − P ( − A ) − ( − A ) (cid:62) P − C = 0 , has a (trivial) solution P = P (cid:62) ∈ R n × n . Now since C (cid:48) = 0 ≥ C , it follows from [42, Theorem 2.3], that theRiccati equation,2 X − X ( − A ) − ( − A ) (cid:62) X =2 X − X ( − A ) − ( − A ) (cid:62) X − C (cid:48) = 0 , has a maximal solution X = P + , which satisfies P + (cid:23) P .So, we have proved that there exists a solution P + to,2 P + P + A + A (cid:62) P + = 0 , which has the property that P + (cid:23) P for any feasiblevariable of (A1). Let ˆ P be the optimal decision variableof (A1). Then, since P + is also feasible for (A1),tr( ˆ P ) ≥ tr( P + ) . (A2)However, P + (cid:23) ˆ P ⇒ ( P + ) ii ≥ ( ˆ P ) ii ⇒ ( P + − ˆ P ) ii ≥ . (A3)By (A2), tr( P + − ˆ P ) ≤ ⇒ n (cid:88) i =1 ( P + − ˆ P ) ii (cid:124) (cid:123)(cid:122) (cid:125) all ≥ ≤ ⇒ ( P + ) ii = ( ˆ P ) ii , i = 1 , ..., n. Hence, P + (cid:23) ˆ P and both have the same diagonal entries.Therefore P + = ˆ P .Consequently, the solution ˆ P to (A1) is the maximal so-lution to ˆ P A + A (cid:62) ˆ P + 2 ˆ P = 0 (cid:3) Remark 1
In practice, for linear systems our fullmethod (15) actually implements,maximize
P,c , [ (cid:15) i ] n (cid:88) i =1 (cid:15) i (A4a) subject to x (cid:62) P x + c ≥ n (cid:88) i =1 (cid:15) i x i , (A4b) P ( A + P ) (cid:22) . (A4c) While (A4) this is not entirely equivalent to (A1) , (A4b) ⇒ (cid:80) i (cid:15) i ≤ tr ( P ) . The implemented method thereforeprovides a very close approximation that also generalizesto the nonlinear case. Appendix B: Iterative improvement
In section III D we displayed the iterative optimizationapproach through which an improved Lyapunov functioncould be obtained, given a suitable initial guess. Theoptimization procedure may be justified by proving thefollowing lemma.
Lemma 2
Let U : R n → R be fixed. Suppose that thereexists α ≥ and V : R n → R such that V ≥ , (B1a) M V (cid:23) , (B1b) ∇ V · ( f + 2 ∇ U ) ≥ α ( f · ∇ U ) + (1 + α ) ||∇ U || (B1c) Then α ( ∇ U · f U ) + ||∇ V − ∇ U || ≤ ∇ V · f V ≤ . Proof . First lower bound ∇ V · f V : ∇ V · f V = ∇ V · ( f + ∇ V )= ∇ V · ( f + 2 ∇ U ) + ||∇ V || − ∇ V · ∇ U (by (B1c)) ≥ α ( f · ∇ U + ||∇ U || ) + ||∇ U || + ||∇ V || − ∇ V · ∇ U = α ( ∇ U · f U ) + ||∇ U − ∇ V || . That ∇ V · f V ≤ (cid:3) Remark 2 (i) If U satisfies the standard Lyapunov con-ditions U ≥ , M U (cid:23) then the choice ( α = 1 , V = U )also satisfies conditions (B1) . Hence, there is always afeasible point which satisfies these inequalities.(ii) If ( α , V ) satisfiers the constraints with < α < ,then α ( ∇ U · f U ) + ||∇ V − ∇ U || ≤ ∇ V · f V ≤ indicates that an improved lower-bound has been foundfor V · f V . The improvement over U · f U is determinedby the difference ||∇ V − ∇ U || and the scaling factor α . (iii) If α = 0 and V satisfies the constraints, then (B1c) reads ||∇ U − ∇ V || ≤ ∇ V · f V ≤ which forces V = U and indicates that U · f U = 0 , i.e.the original Lyapunov function was a good choice. [1] L. Gammaitoni, P. H¨anggi, P. Jung, and F. Marchesoni,“Stochastic resonance,” Rev. Mod. Phys. , 223–287(1998).[2] C. Gardiner, Stochastic Methods: A Handbook For TheNatural And Social Sciences (Springer, 2009).[3] G. Rigas, A. S. Morgans, R. D. Brackston, and J. F.Morrison, “Diffusive dynamics and stochastic models ofturbulent axisymmetric wakes,” J. Fluid Mech. , R2(2015).[4] R. D. Brackston, J. M. Garc`ıa de la Cruz, A. Wynn,G. Rigas, and J. F. Morrison, “Stochastic modelling andfeedback control of bistability in a turbulent bluff bodywake,” J. Fluid Mech. , 726–749 (2016).[5] P. G. Wolynes, J. N. Onuchic, and D. Thirumalai,“Navigating the folding routes,” Science , 1619–1620(1995).[6] C. H. Waddington,
The strategy of the genes: a discus-sion of some aspects of theoretical biology. (Allen & Un-win, London, 1957).[7] N. Moris, C. Pina, and A. M. Arias, “Transition statesand cell fate decisions in epigenetic landscapes,” NatureRev. Gen. , 693–703 (2016).[8] M. I. Freidlin, A. D. Wentzell, and J. Sz¨ucs, Random per-turbations of dynamical systems , 3rd ed. (Springer, Hei-delberg [u.a], 2012).[9] J. Wang, K. Zhang, and E. Wang, “Kinetic paths, timescale, and underlying landscapes: A path integral frame-work to study global natures of nonequilibrium systemsand networks,” J. Chem. Phys. , 125103 (2010).[10] M. Heymann and E. Vanden-Eijnden, “The geometricminimum action method: A least action principle on thespace of curves,” Commun. Pure Appl. Math. , 1052–1117 (2008).[11] D. Liu, “A numerical scheme for optimal transition pathsof stochastic chemical kinetic systems,” J. Comput. Phys. , 8672–8684 (2008).[12] R. de la Cruz, R. Perez-Carrasco, P. Guerrero, T. Alar-con, and K. M. Page, “Minimum Action Path TheoryReveals the Details of Stochastic Transitions Out of Os-cillatory States,” Phys. Rev. Lett. , 128102 (2018).[13] Z. Chen and M. Freidlin, “Smoluchowskikramers approx-imation and exit problems,” Stoch. Dyn. , 569–585(2005).[14] M. K. Cameron, “Finding the quasipotential for nongra-dient SDEs,” Phys. D: Nonlin. Phenom. , 1532–1550(2012).[15] D. Dahiya and M. K. Cameron, “Ordered Line Inte-gral Methods for Computing the Quasi-Potential,” J. Sci.Comput. , 1351–1384 (2017).[16] M. K. Cameron, “Construction of the quasi-potential forlinear SDEs using false quasi-potentials and a geometricrecursion,” arXiv:1801.00327 . [17] J. Wang, “Landscape and flux theory of non-equilibriumdynamical systems with application to biology,” Adv.Phys. , 1–137 (2015).[18] P. Zhou and T. Li, “Construction of the landscapefor multi-stable systems: Potential landscape, quasi-potential, A-type integral and beyond,” J. Chem. Phys. , 094109 (2016).[19] J. Wang, L. Xu, and E. Wang, “Potential landscape andflux framework of nonequilibrium networks: Robustness,dissipation, and coherence of biochemical oscillations,”Proc. Nat. Acad. Sci. U.S.A. , 12271–12276 (2008).[20] C. Lv, X. Li, F. Li, and T. Li, “Constructing the En-ergy Landscape for Genetic Switching System Driven byIntrinsic Noise,” PLoS ONE , e88167 (2014).[21] H. Ge, H. Qian, and X. S. Xie, “Stochastic PhenotypeTransition of a Single Cell in an Intermediate Regionof Gene State Switching,” Phys. Rev. Lett. , 078101(2015).[22] C. Li, T. Hong, and Q. Nie, “Quantifying the land-scape and kinetic paths for epithelialmesenchymal tran-sition from a core circuit,” Phys. Chem. Chem. Phys. ,17949–17956 (2016).[23] J. Wang, L. Xu, E. Wang, and S. Huang, “The potentiallandscape of genetic circuits imposes the arrow of time instem cell differentiation,” Biophys. J. , 29–39 (2010).[24] C. Li, E. Wang, and J. Wang, “Potential Landscape andProbabilistic Flux of a Predator Prey Network,” PLoSONE , e17888 (2011).[25] J. Guo, F. Lin, X. Zhang, V. Tanavde, and J. Zheng,“NetLand: quantitative modeling and visualization ofWaddingtons epigenetic landscape using probabilistic po-tential,” Bioinf. , 1583–1585 (2017).[26] J. X. Zhou, M. D. S. Aliyu, E. Aurell, and S. Huang,“Quasi-potential landscape in complex multi-stable sys-tems,” J. R. Soc. Interface , 3539–3553 (2012).[27] J. Jost, Dynamical Systems (Springer, 2005).[28] P. A. Parrilo,
Structured Semidefinite Programs andSemialgebraic Geometry Methods in Robustness and Op-timization , PhD thesis, Caltech (2000).[29] Note that because a is a fixed point, such a path willtake infinite time. However given any small perturbationaway from a , the time will be finite.[30] I. Dunning, J. Huchette, and M. Lubin, “JuMP: A Mod-eling Language for Mathematical Optimization,” SIAMRev. , 295–320 (2017).[31] A. Papachristodoulou, J. Anderson, G. Valmorbida,S. Prajna, P. Seiler, and P. A. Parrilo, SOSTOOLS:Sum of squares optimization toolbox for MATLAB , http://arxiv.org/abs/1310.4716 (2013).[32] D. Silk, P. D. W. Kirk, C. P. Barnes, T. Toni, A. Rose,S. Moon, M. J. Dallman, and M. P. H. Stumpf, “De-signing attractive models via automated identification of chaotic and oscillatory dynamical regimes.” Nature Com-mun. , 489 (2011).[33] K. G. Murty and S. N. Kabadi, “Some NP-completeproblems in quadratic and nonlinear programming,”Math. Programm. , 117–129 (1987).[34] In practice, these inequalities are imposed by specifyingthat U ( x ) − (cid:15) (cid:80) i x i ∈ Σ and −∇ U ( x ) · f ( x ) ∈ Σ, where Σdenotes the set of polynomials that are Sum of Squares.[35] F. Zhang,
The Schur Complement and Its Applications ,Numerical Methods and Algorithms (Springer, Boston,MA, 2005).[36] A. Papachristodoulou and S. Prajna, “On the construc-tion of Lyapunov functions using the sum of squares de-composition,” in , Vol. 3 (2002) pp. 3482–3487.[37] A. Majumdar, A. A. Ahmadi, and R. Tedrake, “Controland verification of high-dimensional systems with DSOS and SDSOS programming,” in (2014) pp. 394–401.[38] For a linear system in which the matrix A is normal, aminimal basis will also be sufficient since sparsity of thematrix is maintained through (27).[39] C. Kwon, P. Ao, and D. J. Thouless, “Structure ofstochastic dynamics near fixed points,” Proc. Nat. Acad.Sci. U.S.A. , 13029–13033 (2005).[40] R. S. Maier and D. L. Stein, “A scaling theory of bifur-cations in the symmetric weak-noise escape problem,” J.Stat. Phys. , 291–357 (1996).[41] This may be verified by checking that ∂f ∂x = ∂f ∂x .[42] I. Gohberg, P. Lancaster, and L. Rodman, “On Hermi-tian Solutions of the Symmetric Algebraic Riccati Equa-tion,” SIAM J. Control Optim.24