The construction of effective and informative landscapes for stochastic dynamical systems has proven a long-standing and complex problem, including for biological systems. Such landscapes may refer to a true energy function for cases such as protein folding or to a phenomenological metaphor in the case of Waddington’s epigenetic landscape. In many situations, constructing a landscape comes down to obtaining the quasi-potential, a scalar function that quantifies the likelihood of reaching each point in the state-space.
In this work we provide a novel method for constructing such landscapes using a tool from control theory: the Sum-of-Squares method for generating Lyapunov functions. Applicable to any system described by polynomials, this method provides an analytical polynomial expression for the potential landscape, in which the coefficients of the polynomial are obtained via a convex optimisation. The resulting landscapes are based upon a decomposition of the vector-field of the original system, such that the inner product of the gradient of the potential and the remaining dynamics is everywhere negative. By satisfying this condition, our derived landscapes provide both upper and lower bounds on the true quasi-potential; these bounds becoming tight if the decomposition is orthogonal. The method is implemented in the programming language Julia, and is demonstrated to correctly compute the quasi-potential for high-dimensional linear systems and also for a number of nonlinear examples. For a multi-stable stochastic system analogous to a developing stem cell, we use the computed potential to evaluate bounds on the relative likelihood of reaching fixed points equivalent to each of the differentiated phenotypes.