Climbing Escher's stairs: a way to approximate stability landscapes in multidimensional systems

TL;DR

This paper introduces a numerical method and an R package to approximate stability landscapes in multidimensional systems by decomposing differential equations into gradient and non-gradient parts, enabling potential landscape estimation where it traditionally does not exist.

Contribution

The authors present a novel numerical approach and software tool to approximate stability landscapes in complex systems lacking scalar potentials.

Findings

The method effectively decomposes differential equations into gradient and non-gradient components.

Approximate potential landscapes can be constructed in regions where non-gradient effects are small.

An R package implementation makes the approach accessible for practical use.

Abstract

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, for most biological 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. Additionally, we introduce a numerical method for decomposing differential equations into two terms: the gradient term that has an associated potential, and the 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.

Figures

Equations

$$\{\begin{array}{l}\frac{dx}{dt}=f(x,y)=-\frac{\partial V(x,y)}{\partial x}\\ \frac{dy}{dt}=g(x,y)=-\frac{\partial V(x,y)}{\partial y}\end{array}$$

$$\frac{\partial f}{\partial y}=\frac{\partial g}{\partial x}$$

$$V(x,y)=V(x_0,y_0)-{\int }_{x_0}^{x}f(\xi ,y_0)d\xi -{\int }_{y_0}^{y}g(x,\eta )d\eta$$

$$\left[\begin{array}{c}f(x,y)\\ g(x,y)\end{array}\right]=\vec{f}(\vec{x})$$

$$\vec{f}(\vec{x})={\vec{f}}_g(\vec{x})+{\vec{f}}_{ng}(\vec{x})$$

$$\vec{f}(\vec{x})\approx \vec{f}(\overrightarrow{x_0})+J(\overrightarrow{x_0})\Delta \vec{x}$$

$$J=J^T$$

$$\{\begin{array}{l}{M}_{symm}=\frac{1}{2}\left(M+M^T\right)\\ M_{skew}=\frac{1}{2}\left(M-M^T\right)\end{array}$$

$$\vec{f}(\vec{x})\approx \vec{f}(\overrightarrow{x_0})+J_{symm}(\overrightarrow{x_0})\Delta \vec{x}+J_{skew}(\overrightarrow{x_0})\Delta \vec{x}$$

$$\{\begin{array}{l}{\vec{f}}_g(\vec{x})\approx \vec{f}(\overrightarrow{x_0})+J_{symm}(\overrightarrow{x_0})\Delta \vec{x}\\ {\vec{f}}_{ng}(\vec{x})\approx J_{skew}(\overrightarrow{x_0})\Delta \vec{x}\end{array}$$

$$\Delta V(\overrightarrow{x_1},\overrightarrow{x_0})\equiv V(\overrightarrow{x_1})-V(\overrightarrow{x_0})\approx -\vec{f}({\vec{x}}_0)\cdot \Delta \vec{x}-\frac{1}{2}\Delta {\vec{x}}^T{J}_{symm}(\overrightarrow{x_0})\Delta \vec{x}$$

$$V({\vec{x}}_{next})\approx V({\vec{x}}_{previous})+\Delta V({\vec{x}}_{next},{\vec{x}}_{previous})$$

$$V(x_i,y_j)=V(x_0,y_0)+\sum _{k=1}^i\Delta V(x_k,y_0;x_{k-1},y_0)+\sum _{l=1}^j\Delta V(x_i,y_l;x_i,y_{l-1})$$

$$err(\vec{x})\approx \frac{\mid {\vec{f}}_{ng}(\vec{x})\mid }{\mid {\vec{f}}_g(\vec{x})\mid +\mid {\vec{f}}_{ng}(\vec{x})\mid }\approx \frac{\mid J_{skew}(\vec{x})\mid }{\mid J_{symm}(\vec{x})\mid +\mid J_{skew}(\vec{x})\mid }$$

$$\{\begin{array}{l}\frac{dx}{dt}=f(x,y)=-x(x^2-1)+p_x(x,y)\\ \frac{dy}{dt}=g(x,y)=-y(y^2-1)+p_y(x,y)\end{array}$$

$$\{\begin{array}{l}\frac{dx}{dt}=-y\\ \frac{dy}{dt}=x\end{array}$$

$$\{\begin{array}{l}\frac{dx}{dt}=b_x-r_{x}x+\frac{a_x}{k_x+y^n}\\ \frac{dy}{dt}=b_y-r_{y}y+\frac{a_y}{k_y+x^n}\end{array}$$

$$\{\begin{array}{l}\frac{dx}{dt}=ax-bxy\\ \frac{dy}{dt}=cxy-dy\end{array}$$

$$\frac{\partial f_i}{\partial x_j}=\frac{\partial f_j}{\partial x_i}:i\ne j$$

$$\frac{d{x}_i}{dt}=f_i(\vec{x})=-\frac{\partial V}{\partial x_i}:i=\mathrm{1..}n$$

$$V(\vec{x})=V({\vec{x}}_0)-{\int }_{\Gamma }\sum _{i=1}^n{f}_i(\vec{x})d{x}_i$$

$$V(x_1,y_0)\approx V(x_0,y_0)+\Delta V(x_1,y_0;x_0,y_0)$$

$$V(x_3,y_0)\approx V(x_0,y_0)+\Delta V(x_1,y_0;x_0,y_0)+\Delta V(x_2,y_0;x_1,y_0)+\Delta V(x_3,y_0;x_2,y_0)$$

$$\begin{array}{r}V(x_3,y_2)\approx V(x_0,y_0)+\Delta V(x_1,y_0;x_0,y_0)+\Delta V(x_2,y_0;x_1,y_0)+\Delta V(x_3,y_0;x_2,y_0)+\\ +\Delta V(x_3,y_1;x_3,y_0)+\Delta V(x_3,y_2;x_3,y_1)\end{array}$$

Full text

1 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).

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.

1.1 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)).

[TABLE]

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)).

[TABLE]

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.

[TABLE]

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)).

[TABLE]

2 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)).

[TABLE]

The gradient term ($\vec{f}_{g}(\vec{x})$) captures the part of the system that can be associated to a potential function, while the non-gradient term ($\vec{f}_{ng}(\vec{x})$) represents the deviation from this ideal case. We’ll use $\vec{f}_{g}(\vec{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 $\mid\vec{f}_{ng}(\vec{x})\mid$. In regions where the gradient term is stronger than the non-gradient 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 $\vec{f}(\vec{x})$ can be approximated around a point $\vec{x_{0}}$ using equation (6), where $J(\vec{x_{0}})$ is the Jacobian matrix evaluated at the point $\vec{x_{0}}$ and $\Delta\vec{x}$ is defined as the distance to this point, that is, $\Delta\vec{x}=\vec{x}-\vec{x_{0}}$, written as a column vector.

[TABLE]

As usual in linearization, we have neglected the terms of order $2$ and higher in equation (6). This approximation is valid for $\vec{x}$ close to $\vec{x_{0}}$. For the system to be gradient (equation (2)) its Jacobian has to be symmetric (see equation (7), where $T$ represents transposition).

[TABLE]

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)).

[TABLE]

Using the skew-symmetric decomposition described in equation (8), we can rewrite (6) as:

[TABLE]

Equation (9) represents a natural, well-defined and operational way of writing our vector field $\vec{f}(\vec{x})$ decomposed as in equation (5), that is, as the sum of a gradient and a non-gradient term (see (10)).

[TABLE]

The gradient term $f_{g}(\vec{x})$ can thus be associated to a potential $V(\vec{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 $\vec{x_{0}}$ and another point $\vec{x_{1}}\equiv\vec{x_{0}}+\Delta\vec{x}$ separated by a distance $\Delta\vec{x}$ (see equation (11)).

[TABLE]

Provided we know the value of the potential at one point $\vec{x}_{previous}$, equation (11) allows us to estimate the potential at a different point $\vec{x}_{next}$ (cf.: equation (12)).

[TABLE]

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.

[TABLE]

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:

It has been derived from a set of linearized equations, sampled over a grid 2. 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 $\mid\Delta\vec{x}\mid^{3}$, where $\mid\Delta\vec{x}\mid$ 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 $\vec{f}_{ng}(\vec{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).

[TABLE]

2.1 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).

3 Results

3.1 Synthetic examples

3.1.1 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)).

[TABLE]

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\cdot m(x,y)$ and $p_{y}(x,y)=-0.4x\cdot 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.

3.1.2 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.

[TABLE]

3.2 Biological examples

3.2.1 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)).

[TABLE]

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$, $foldXY=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).

3.2.2 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.

[TABLE]

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).

4 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.

5 Acknowledgments

This work was greatly inspired by the discussions with Cristina Sargent, Iñaki Úcar, Enrique Benito, Sanne J.P. van den Berg, Tobias Oertel-Jäger, Jelle Lever, and Els Weinans. This work was supported by funding from the European Union’s Horizon 2020 research and innovation programme for the ITN CRITICS under Grant Agreement Number 643073.

Online appendix

A.1 Gradient conditions for a system with an arbitrary number

of dimensions

Dynamics in equation (1) and the condition for the crossed derivatives (2) can be straightforwardly generalized (see equations (A.2) and (A.1) to systems with an arbitrary number of state variables $\vec{x}=(x_{1},...,x_{n})$. Particularly, if and only if our system of equations $\frac{dx_{i}}{dt}=f_{i}(\vec{x})$ satisfies the condition:

[TABLE]

then a potential $V(\vec{x})$ exists related to the original vector field:

[TABLE]

and such a potential can be computed using a line integral:

[TABLE]

where the line integral in (A.3) is computed along any curve $\Gamma$ joining the points $\vec{x}_{0}$ and $\vec{x}$.

It is important to note that the number of equations ($N$) described in condition (A.1) grows rapidly with the dimensionality of the system ($D$), following the series of triangular numbers $N=\frac{1}{2}(D-1)D$. Thus, the higher the dimensionality, the harder it may get to fulfill condition (A.1). As a side effect, we see that one-dimensional systems have zero conditions and their stability landscape is thus always well-defined.

A.2 Detailed example of application

To calculate the value of $V$ at, for instance, the point $(x_{3},y_{2})$ of a grid, we should begin by assigning [math] to the potential at our arbitrary starting point (i.e.: $V(x_{0},y_{0})=0$ by definition). Then, we need a trajectory that goes from $(x_{0},y_{0})$ to $(x_{3},y_{2})$, iterating over the intermediate grid points (see figure A.1).

In the first step we go from $(x_{0},y_{0})$ to $(x_{1},y_{0})$. The new potential is thus (using (12)):

[TABLE]

The next two steps continue in the horizontal direction, all the way to $(x_{3},y_{0})$. The value of the potential there is:

[TABLE]

Now, to reach our destination $(x_{3},y_{2})$ we have to move two steps in the vertical direction:

[TABLE]

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.