Climbing Escher’s stairs: A way to approximate stability landscapes in multidimensional systems

Stability landscapes are useful for understanding the properties of dynamical systems. These landscapes can be calculated from the system’s dynamical equations using the physical concept of scalar potential. Unfortunately, it is well known that for most systems with two or more state variables such potentials do not exist. Here we use an analogy with art to provide an accessible explanation of why this happens and briefly review some of the possible alternatives. Additionally, we introduce a novel and simple computational tool that implements one of those solutions: the decomposition of the differential equations into a gradient term, that has an associated potential, and a non-gradient term, that lacks it. In regions of the state space where the magnitude of the non-gradient term is small compared to the gradient part, we use the gradient term to approximate the potential as quasi-potential. The non-gradient to gradient ratio can be used to estimate the local error introduced by our approximation. Both the algorithm and a ready-to-use implementation in the form of an R package are provided.


Introduction
With knowledge becoming progressively more interdisciplinary, the relevance of science communication is rapidly increasing.Mathematical concepts are among the hardest topics to communicate to non-expert audiences, policy makers, and also to scientists with little mathematical background.Visual methods are known to be successful ways of explaining mathematical concepts and results to non-specialists.
One particularly successful visualization method is that of the stability landscape, also known as the rolling marble or ball-in-a-cup diagram Edelstein-Keshet (1988), Strogatz (1994), Beisner, Haydon, and Cuddington (2003), Pawlowski (2006), which origin can be traced back to the introduction of the scalar potential in physics by Lagrange in the 18th century Lagrange (1777).In stability landscapes (e.g.: figure 1) the horizontal position of the marble represents the state of the system at a given time.With this picture in mind, the shape of the surface represents the underlying dynamical rules, where the slope is proportional to the speed of the movement.The peaks on the undulated surface represent unstable equilibrium states and the wells represent stable equilibria.Different basins of attraction are thus separated by peaks in the surface.Stability landscapes have proven to be a successful tool to explain advanced concepts related with the stability of dynamical systems in an intuitive way.Some examples of those advanced concepts are multistability, basin of attraction, bifurcation points and hysteresis (see Scheffer et al. (2001), Beisner, Haydon, and Cuddington (2003) and figure 1).
The main reason for the success of this picture arises from the fact that stability landscapes are built as an analogy with our most familiar dynamical system: movement of a marble.Particularly, the movement of a marble along a curved landscape under the influence of its own weight.The stability landscape corresponds then with the physical concept of potential energy Strogatz (1994).This explains why our intuition, based in what we know about movement in our everyday life, works so well reading this diagrams.It is important to stress the fact that under this picture there's not such a thing as inertia Pawlowski (2006).The accurate analogy is that of a marble rolling in a surface while submerged inside a very viscous fluid Strogatz (1994).Like with any other analogy, it is important to be aware of its limitations.The one we address here is the fact that, for models with more than one state variable, such a potential doesn't exist in general.To get an intuitive feeling of why this is true, picture a model with a stable cyclic attractor.As the slope of the potential should reflect the speed of change, we would need a potential landscape where our marble can roll in a closed loop while always going downhill.Such a surface is a classical example of an impossible object (see figure 2 and L. S. Penrose and Penrose (1958) for details).
Figure 2: The Penrose stair L. S. Penrose and Penrose (1958) is a classical example of an impossible object.In such a surface, it is possible to walk in a closed loop while permanently going downhill.The scalar potential of a system with a cyclic attractor, if existed, should have the same impossible geometry.This object was popularized by the Dutch artist M.C.Escher (for two beautiful examples, see Escher (1960) and Escher (1961)).
As this is a centuries-old problem (see for instance Helmholtz (1858)), it is perhaps of not surprise that several methods have been proposed to approximate stability landscapes for general, high-dimensional systems.For a comprehensive mathematical review, we refer to Zhou et al. (2012).Here, different criteria are used to decompose analytically the dynamics in a scalar potential and an extra term.Some interesting alternatives are presented in Pawlowski (2006), like potentials for second-order systems or the use of Lyapunov functions as stability landscapes.For stochastic systems, the Freidlin-Wentzell potential Freidlin and Wentzell (2012) has been proposed as a very strong candidate due to its links with transition rates Zhou et al. (2012), Nolting and Abbott (2015).
In the present work we review the state of the art, and we introduce a simple new method to deal with the fundamental problem of approximating stability landscapes for high dimensional systems.Specificially, we introduce an algorithm to decompose differential equations as the sum of a gradient and a non-gradient part.Each part can be used, respectively, to compute an associated potential and to measure the local error introduced by our picture.In order to reach those interested readers with little background in mathematics, we limited our mathematical weaponry.Knowledge of basic linear algebra and basic vector calculus will suffice to understand the paper to its last detail.Additionally, we provide a ready to use R package that implements the algorithm this paper describes.

Mathematical background
Consider a coupled differential equation with two state variables x and y (see online appendix A.1 for a generalization in more dimensions).The dynamics of such system can be described by a stability landscape if a single two-dimensional function V (x, y) exists whose slope is proportional to the change in time of both states (equation ( 1)).
It can be shown that such potential V (x, y) only exists if the crossed derivatives of functions f (x, y) and g(x, y) are equal for all x and y (equation ( 2)).Systems satisfying equation (2) are known as conservative, irrotational or gradient fields (cf.section 8.3 of Marsden and Tromba (2003)).
If condition (2) holds, we can use a line integral (Marsden and Tromba (2003), section 7.2) to invert (1) and calculate V (x, y) using the functions f (x, y) and g(x, y) as an input.An example of this inversion is equation ( 3), where we have chosen an integration path composed of a horizontal and a vertical line.
The attentive reader may have raised her or his eyebrow after reading the word chosen applied to an algorithm.In fact, we can introduce this arbitrary choice without affecting the final result.The condition for potentials to exist (equation ( 2)) implies that any line integral between two points in this vector field should be independent of the path (cf.section 7.2 of Marsden and Tromba (2003)).Going back to our rolling marble analogy, we can gain some intuition about why this is true: in a landscape the difference in potential energy between two points is proportional to the difference in height, and thus stay the same for any path.If the condition (2)) is not fulfilled, the potential calculated with (3)) will depend on the chosen integration path.As this is an arbitrary choice, the computed potential will be an artifact with no natural meaning.
For the sake of a compact and easy to generalize notation, in the rest of this paper we will arrange the equations of our system as a column vector (equation ( 4)).

Methods
The method for deriving a potential we propose is based on the decomposition of a vector field in a conservative or gradient part and a non-gradient part (see equation ( 5)).
The gradient term ( f g ( x)) captures the part of the system that can be associated to a potential function, while the non-gradient term ( f ng ( x)) represents the deviation from this ideal case.We'll use f g ( x) to compute an approximate or quasi-potential.The absolute error of this approach is given as the euclidean size of the non-gradient term | f ng ( x) |.In regions where the gradient term is stronger than the nongradient term, the condition (2) will be approximately fulfilled, and thus the calculated quasi-potential will represent an acceptable approximation of the underlying dynamics.Otherwise, the non-gradient term is too dominant to approximate a potential landscape.
In order to achieve a decomposition like (5), we begin by linearizing our model equations.Any sufficiently smooth and continuous vector field f ( x) can be approximated around a point x 0 using equation ( 6), where J( x 0 ) is the Jacobian matrix evaluated at the point x 0 and ∆ x is defined as the distance to this point, that is, ∆ x = x − x 0 , written as a column vector.
As usual in linearization, we have neglected the terms of order 2 and higher in equation ( 6).This approximation is valid for x close to x 0 .For the system to be gradient (equation ( 2)) its Jacobian has to be symmetric (see equation ( 7), where T represents transposition).
We know from basic linear algebra that any square matrix M can be uniquely decomposed as the sum of a skew and a symmetric matrix (see equation ( 8)).
Using the skew-symmetric decomposition described in equation ( 8), we can rewrite (6) as: Equation ( 9) represents a natural, well-defined and operational way of writing our vector field f ( x) decomposed as in equation ( 5), that is, as the sum of a gradient and a non-gradient term (see ( 10)).
The gradient term f g ( x) can thus be associated to a potential V ( x) that can be computed analytically for this linearized model using a line integral (see equation (3) for the two dimensional case, or (A.3) in the online appendix for the general one).The result of this integration yields an analytical expression for the potential difference between the reference point x 0 and another point x 1 ≡ x 0 + ∆ x separated by a distance ∆ x (see equation ( 11)).
∆V ( Provided we know the value of the potential at one point x previous , equation ( 11) allows us to estimate the potential at a different point x next (cf.: equation ( 12)).
Equation ( 12) can be applied sequentially over a grid of points to calculate the approximate potential on each of them.In two dimensions, the resulting potential is given by the closed formula (13).The cases with 3 and more dimensions can be generalized straightforwardly.For a step by step example, see section A.2 in the online appendix.For a flowchart overview of the algorithm, please refer to figure 3.
As with any other approximation we need a way to estimate and control its error.The stability landscape described in (11) has two main sources of errors: 1.It has been derived from a set of linearized equations, sampled over a grid 2. It completely neglects the effects of the non-gradient part of the system The error due to linearization in equation ( 11) is roughly proportional to | ∆ x | 3 , where | ∆ x | is the euclidean distance to the reference point.By introducing a grid, we expect the linearization error to decrease with the grid's step size.
The more fundamental error due to neglecting the non-gradient component of our system is not affected by the grid's step choice.From equation ( 5) it is apparent that we can use the euclidean size of f ng ( x) as an approximation of the local error introduced by our algorithm.The relative error due to this effect can be estimated using equation ( 14).

Implementation
As an application of the above-mentioned ideas, and following the spirit of reproducible research, we developed a transparent R package we called rolldown Rodríguez-Sánchez (2019).Our algorithm accepts a set of dynamical equations and a grid of points defining our region of interest as an input.The output is the estimated potential and the estimated error, both of them calculated at each point of our grid (see figure 3 for details).

A four well potential
We first tested our algorithm with a synthetic model of two uncoupled state variables.Uncoupled systems are always gradient as all non-diagonal values of the Jacobian are zero everywhere.We added the interaction terms p x and p y to be able to make it gradually non-gradient (see equation ( 15)).
When we choose those non-gradient interactions to be zero, the system is purely gradient and corresponds with a four-well potential.Our algorithm rendered it successfully and with zero error (cf.figure 4, row A).
In order to test our algorithm, we introduced a non-gradient interaction of the form p x (x, y) = 0.3y •m(x, y) and p y (x, y) = −0.4x• m(x, y), with m(x, y) = e (x−1) 2 +(y−1) 2 .m(x, y) serves as a masking function guaranteeing that our interaction term will be negligible everywhere but in the vicinity of (x, y) = (1, 1).
After introducing this non-gradient interaction a four-well potential is still recognizable (cf.figure 4, row B).As expected, the error was correctly captured to be zero everywhere but in the region around (x, y) = (1, 1).The error map warns us against trusting the quasi-potential in the upper right region, and guarantees that elsewhere it will work fine.Notice that, accordingly, the upper right stable equilibrium falls slightly outside its corresponding well.The rest of equilibria, to the contrary, fit correctly inside their corresponding wells.Row A shows the application to a gradient case (equation ( 15) with interaction terms equal to zero).As expected, the error is zero everywhere.In row B our algorithm is applied to a non-gradient case (equation ( 15), with non-zero interaction terms).

A fully non-gradient system
In order to stress our algorithm to the maximum, we tested it in a worst case scenario: a system with zero gradient part everywhere.Particularly, we fed it with equation ( 16).All solutions (but the unstable equilibrium point at (0, 0)) are cyclic (cf.5, left panel)).As we discussed in the introduction, calculating a potential for a system with cyclic trajectories is as impossible as Escher's paintings (and for similar reasons).This fact is captured by our algorithm, that correctly predicts a relative error of 1 everywhere (see figure 5, central panel).In this case, the quasi-potential (figure 5, left panel) is not even locally useful.16)).In all panels the dots represent center equilibrium points.The left panel shows the phase plane containing the actual "deterministic skeleton" of the system.The central panel shows the quasi-potential.The right panel shows the estimated error.

A simple regulatory gene network
Waddington's epigenetic landscapes Gilbert (1991) are a particular application of stability landscapes to gene regulatory networks controlling cellular differentiation.When applied to this problem, stability/epigenetic landscapes serve as a visual metaphor for the branching pathways of cell fate determination Bhattacharya, Zhang, and Andersen (2011).
A bistable network cell fate model can be described by a set of equations like (17).Such a system represents two genes (x and y) that inhibit each other.This circuit works as a toggle switch with two stable steady states, one with dominant x, the other with dominant y (see Bhattacharya, Zhang, and Andersen (2011)).Our parameter choice (a x = 0.4802, a y = 0.109375, b x = 0.2, b y = 0.3, k x = 0.2401, k y = 0.0625, r x = r y = 1 and , n = 4) corresponds with equations 6 and 7 of Bhattacharya, Zhang, and Andersen (2011), where we modified two parameters (By = 0.3, f oldXY = 1.75, in their notation) in order to induce an asymmetry in the the dynamics.Despite this system is non-gradient, our algorithm correctly predicts the existence of two wells (see figure 6, row D).The relative error, despite being distinct from zero in some regions, is not very high.This means that our quasi-potential is a reasonable approximation of the underlying dynamics.Indeed, the equilibria correspond to the wells (stable) and the peak (unstable).

Predator prey dynamics
The Lotka-Volterra equations ( 18) are a classical predator-prey model Volterra (1926).In this model x and y represent, respectively, the prey and predator biomasses.This model is known to have cyclic dynamics.As we discussed in our analogy with Escher's paintings, we cannot compute stability landscapes in the regions of the phase plane where the dynamics are cyclic.When we apply our method to a system like equation ( 18), the error map correctly captures the fact that our estimated potential is not trustworthy (see figure 6, row E).In row D we applied our algorithm to the simple gene regulatory network described in equation ( 17).
In row E we apply our algorithm to a Lotka-Volterra system (equation ( 18), with a = b = c = d = 1).

Discussion
The use of stability landscapes as a helping tool to understand one-dimensional dynamical systems achieved great success, especially in interdisciplinary research communities.A generalization of the idea of scalar potential to two-dimensional systems seemed to be a logical next step.Unfortunately, as we have seen, there are reasons that make two (and higher) dimensional systems fundamentally different from the one-dimensional case.The generalization, straightforward as it may look, is actually impossible for most dynamical equations.A good example of this impossibility is any system with cyclic dynamics, whose scalar potential should be as impossible a the Penrose stairs in Escher's paintings.
As a consequence, any attempt of computing stability landscapes for high-dimensional systems should, necessarily, drop some of the desirable properties of classical scalar potentials.For instance, the method proposed by Bhattacharya Bhattacharya, Zhang, and Andersen (2011) smartly avoids the problem of path dependence of line integrals by integrating along trajectories, removing thus the freedom of path choice.The price paid is that Bhattacharya's algorithm cannot guarantee continuity along separatrices between basins of attraction for general multistable systems.
We share the perception with other authors Pawlowski ( 2006) that the concept of potential is often misunderstood in research communities with a limited mathematical background.Analytical methods like the ones described in Zhou et al. (2012) or Nolting and Abbott (2015) require familiarity with advanced mathematical concepts such as partial differential equations.The algorithm we present here is an attempt to preserve as much as possible from the classical potential theory while keeping the mathematical complexity as low as possible.Additionally, our algorithm provides: • Integrity.At each step the strength of the non-gradient term is calculated.If this term is high, it is fundamentally impossible to calculate a scalar potential with any method.If this term is zero, our solution converges to the classical potential.• Safety.The relative size of the non-gradient term can be interpreted as an error term, mapping which regions of our stability landscape are dangerous to visit.• Speed.The rendering of a printing quality surface can be performed in no more than a few seconds in a personal laptop.• Simplicity.The required mathematical background is covered by any introductory course in linear algebra and vector calculus.• Generality.The core of the algorithm is the skew-symmetric decomposition of the Jacobian.This operation can be easily applied to square matrices of any size, generalizing our algorithm for working in 3 or more dimensions.• Usability.We provide the algorithm in the form of a ready to use, documented and tested R package.
It is important to notice that, despite our algorithm provides us with a way of knowing which regions of the phase plane can be "safely visited", we cannot navigate the phase plane freely but only along trajectories.This interplay between regions and trajectories limits the practical applicability of our algorithm to those trajectories that never enter regions with high error.
The concept of potential is paramount in physical sciences.The main reason for the ubiquity of potentials in physics is that several (idealized) physical systems are known to be governed only by gradient terms (e.g.: movement in friction-less systems, classical gravitatory fields, electrostatic fields, . . .).As physical potentials can be related with measurable concepts like energy, its use goes way further than visualization.
From the depth and width of a potential we can learn about transition rates and resilience to pulse perturbations.The height of the lowest barrier determines the minimum energy to transition to an alternative stable state, which relates to the probability of a critical transition in a stochastic environment Zhou et al. (2012), Nolting and Abbott (2015), Hänggi, Talkner, and Borkovec (1990).All these results remain true for non-physical problems that happen to be governed exclusively by gradient dynamics, and, we claim, should remain approximately true for problems governed by weakly non-gradient dynamics.This is the situation our algorithm has been designed to deal with.
Regarding visualization alone, it may be worth reconsidering why do we prefer the idea of stability landscape over a traditional phase plane figure, especially after pointing out all the difficulties of calculating stability landscapes for higher-dimensional systems.It is true that the phase plane is slightly less intuitive than the stability landscape, but it has a very desirable property: it doesn't require the imagination of a surrealist artist to exist.
In the first step we go from (x 0 , y 0 ) to (x 1 , y 0 ).The new potential is thus (using ( 12)): V (x 1 , y 0 ) ≈ V (x 0 , y 0 ) + ∆V (x 1 , y 0 ; x 0 , y 0 ) The next two steps continue in the horizontal direction, all the way to (x 3 , y 0 ).The value of the potential there is: V (x 3 , y 0 ) ≈ V (x 0 , y 0 ) + ∆V (x 1 , y 0 ; x 0 , y 0 ) + ∆V (x 2 , y 0 ; x 1 , y 0 ) + ∆V (x 3 , y 0 ; x 2 , y 0 ) Now, to reach our destination (x 3 , y 2 ) we have to move two steps in the vertical direction: V (x 3 , y 2 ) ≈ V (x 0 , y 0 ) + ∆V (x 1 , y 0 ; x 0 , y 0 ) + ∆V (x 2 , y 0 ; x 1 , y 0 ) + ∆V (x 3 , y 0 ; x 2 , y 0 )+ +∆V (x 3 , y 1 ; x 3 , y 0 ) + ∆V (x 3 , y 2 ; x 3 , y 1 ) Generalizing the previous example we see that we can compute the approximate potential at a generic point (x i , y j ) using the closed formula (13).Both our example (A.2) and formula ( 13) have been derived sweeping first in the horizontal direction and next in the vertical one.Of course, we can choose different paths of summation.Nevertheless, because we are building our potential neglecting the non-gradient part of our vector field, we know that our results will converge to the same solution regardless of the chosen path.

Figure 1 :
Figure 1: Example of a set of 4 stability landscapes (reproduced from Scheffer et al. (2001)) used to illustrate bistability in ecosystems (e.g: forest/desert, eutrophicated lake/clear water, etc.).The upper side of the figure shows the stability landscape of a one-dimensional system for 4 different values of a control parameter.The lower side shows the bifurcation diagram.This diagram proved to be a successful tool for explaining advanced concepts in dynamical systems theory such as bistability and fold bifurcations to scientific communities as diverse as ecologists, mathematicians and environmental scientists.

Figure 3 :
Figure 3: Flowchart showing the basic functioning of our implementation of the algorithm described in this paper.

Figure 4 :
Figure4: Results for two synthetic examples.In all panels the dots represent equilibrium points (black for stable, otherwise white).The left panel shows the phase plane containing the actual "deterministic skeleton" of the system.The central panel shows the quasi-potential.The right panel shows the estimated error.Row A shows the application to a gradient case (equation (15) with interaction terms equal to zero).As expected, the error is zero everywhere.In row B our algorithm is applied to a non-gradient case (equation (15), with non-zero interaction terms).

Figure 5 :
Figure5: Results for a fully non-gradient system (equation (16)).In all panels the dots represent center equilibrium points.The left panel shows the phase plane containing the actual "deterministic skeleton" of the system.The central panel shows the quasi-potential.The right panel shows the estimated error.
dx dt = b x − r x x + ax kx+y n dy dt = b y − r y y +

Figure 6 :
Figure6: Results for two biological systems.In all panels the dots represent equilibrium points (black for stable, otherwise white).The left panel shows the phase plane containing the actual "deterministic skeleton" of the system.The central panel shows the quasi-potential.The right panel shows the estimated error.In row D we applied our algorithm to the simple gene regulatory network described in equation (17).In row E we apply our algorithm to a Lotka-Volterra system (equation (18), with a = b = c = d = 1).