## Abstract

The developmental dynamics of multicellular organisms is a process that takes place in a multi-stable system in which each attractor state represents a cell type, and attractor transitions correspond to cell differentiation paths. This new understanding has revived the idea of a quasi-potential landscape, first proposed by Waddington as a metaphor. To describe development, one is interested in the ‘relative stabilities’ of *N* attractors (*N* > 2). Existing theories of state transition between local minima on some potential landscape deal with the exit part in the transition between two attractors in pair-attractor systems but do not offer the notion of a global potential function that relates more than two attractors to each other. Several *ad hoc* methods have been used in systems biology to compute a landscape in non-gradient systems, such as gene regulatory networks. Here we present an overview of currently available methods, discuss their limitations and propose a new decomposition of vector fields that permits the computation of a quasi-potential function that is equivalent to the Freidlin–Wentzell potential but is not limited to two attractors. Several examples of decomposition are given, and the significance of such a quasi-potential function is discussed.

## 1. Introduction: the need for a quasi-potential function

For high-dimensional, nonlinear dynamical systems that exhibit a large number of stable steady states (attractors) far from thermodynamic equilibrium, a characteristic system behaviour of interest is manifest at the scale of the transitions between these attractors across large distances in phase space. Such systems are epitomized by gene regulatory networks of *N* genes, where attractor states represent the cell types in the metazoan body. Accordingly, transitions between attractor states correspond to cell phenotype changes during normal development, homeostasis and disease [1,2]. But standard analysis of dynamical systems focuses on the existence and local properties of a given steady state (linear stability). Thus, it fails to capture the essence of complex systems that operate at time scales where the transitions between attractors constitute the characteristic behaviour.

Specifically, in a system with *M* attractor states , one would like to obtain a sense of the ‘relative depth’ [3] or, more precisely, the ordering of the metastable attractor states through some energy-like function *U*_{1}, *U*_{2}, *U*_{3}, … , *U _{M}*. Such a function

*U*(

**x**) for any point

**x**in the state space of the system would inform us about the probability and direction of transitions between attractor states in a noisy or perturbed system—for instance, in understanding spontaneous cell fate choice (from among a set of alternative nearby attractors for a given initial state) or the feasibility of artificial reprogramming between particular cell types. Such ordering would also explain hierarchies and the arrow of time of ontogenesis.

Let us consider a deterministic system (network) of *N* variables *x _{i}* (e.g. the activity of interacting genes) whose values describe the cell state

**(**

*x**t*) = (

*x*

_{1}(

*t*),

*x*

_{2}(

*t*), … ,

*x*(

_{N}*t*))

^{T}(all notations are given in table 1) and whose dynamics results from how each gene influences the activity of other genes (as invariantly predetermined by the gene regulatory network that is hard-coded in the genome). Such dynamics is described by the first-order ordinary differential equations (ODEs), which in general are nonlinear: 1.1or in vector form: Such a system represents an

*influence network*in which the variables (network nodes) influence each other's values according to the network rules defined by the vector

**F**(

**). Thus,**

*x***F**(

**x**) represents the ‘forces’ acting to change the system state

**(**

*x**x*

_{1},

*x*

_{2}, … ,

*x*) in this inertia-free system. In one-dimensional systems, one always obtains a potential function

_{N}*U*

^{int}(

**x**) by integration, which is a measure of the metastability of two (attractor) states

*x*and

_{A}*x*in a system: 1.2

_{B}Equivalently, the driving force is proportional to the gradient of this potential function: 1.3

In the system at thermodynamic equilibrium where the steady-state probability density function over **x** (giving the presence of thermal fluctuations) satisfies the Boltzmann distribution, *U*^{int} is the potential energy and is the path-independent potential difference between the two states *x** _{A}* and

**. It is related to the transition probability via the Arrhenius equation, 1.4here,**

*x*_{B}*ɛ*is the magnitude of noise [4].

Most works refer to transitions between equilibrium states of thermodynamic systems or simple cases such as Markovian one-dimensional non-equilibrium systems. Here we are interested in the rates of the transition between states in far-from-equilibrium, non-conservative and high-dimensional open systems with a driving force emanating from the internal interactions as described by equation (1.1). For such cases, even if we have a gradient system where a potential function *U* can be obtained by integration, i.e. where there exists a function *U*^{int}(** x**) with the following properties:
1.5such that one can obtain, by integration, a potential function

*U*

^{int}(

**) 1.6even then the transition rate for**

*x*

*x**→*

_{A}

*x**is not path-independent and is determined by the potentials of attractor states and saddle points between them. The transition rate is related to and the transition follows the*

_{B}*least action path*(LAP) [5] as discussed below. Here

*x*

_{s}is the saddle point between two attractors

**x**

*and*

_{A}**x**

*, as shown in figure 1.*

_{B}Yet, a potential function in non-gradient systems that is loosely referred to as ‘quasi-potential’ (indicated here by *Ũ* when its computation is not further specified) is not without significance if properly defined such that it can serve as a quantity for the metastability ordering of (metastable) attractor states of the system in equation (1.1). Specifically, for *Ũ* to be of meaning for our purposes in characterizing cell state transitions in development, we require that *Ũ*

(i) satisfies d

*Ũ*/d*t*< 0 for≠*x** and d*x**Ũ*/d*t*= 0 for=*x**, where*x** are the stable steady states, which can be either fixed points or a limit cycle, according to the stability theory of Lyapunov [6];*x*(ii) is related to the Freidlin–Wentzell action along the LAP between two states and hence permits the computation of the transition rate between them [5];

(iii) we also propose that the following transitivity condition is satisfied: if

*P*_{a}_{→b}>*P*_{b}_{→a}and*P*_{b}_{→c}>*P*_{c}_{→b}then*P*_{a}_{→c}>*P*_{c}_{→a}for any three states*x*,_{a}*x*,_{b}*x*. Here_{c}*P*_{i}_{→j}is the transition probability from a state*x*to another state_{i}*x*._{j}

Note that (i) is the stability property of an attractor state in a non-local sense; (ii) refers to the relationship of any two points; and only (iii) pertains to multiple attractor states. This last requirement (iii) is only conjectured to be satisfied in a subset of systems and will be elaborated in §6. This last condition must be satisfied for *Ũ* to permit global ordering of metastability as described above.

## 2. The physics problem: the quasi-potential function

Because high-dimensional, non-equilibrium systems generally are not gradient systems (i.e. equation (1.5) is not satisfied), 2.1

By contrast, one can enforce a partial notion of a quasi-potential *Ũ*, if we write the driving force as a sum of two terms,
2.2Where **F**_{r} is the ‘remainder’ beyond the component of the driving force that follows the potential gradient. Thus, we decompose the non-gradient vector field **F**(**x**), which is finite and smooth (at least twice differentiable), into two components: one that is the gradient of some ‘potential’ *Ũ* and a second that represents the remainder of driving forces of the dynamic system. The question then is: what is the physical meaning of these terms? Given that there are infinite ways of decomposing a vector field into a sum of two fields, uniqueness of decomposition is imparted by the introduction of constraints, which in turn embody the physical significance. Our objective here is to find decomposition such that the quasi-potential *Ũ* presents exactly the transition barriers associated with the state transitions among steady states and the ‘remainder’ of the driving force **F**_{r} will not contribute to the efforts needed for the transition.

Little work has been done to the construction of quasi-potentials given a system d*x*/d*t* = **F**(**x**) to determine the metastability ordering as defined above. Prigogine [7] proposed a measure of global stability in nonlinear dynamic systems far from equilibrium to compare transition barriers between any pairs of attractors in a multi-stable system following a larger perturbation. Efimov [8] studied the mathematical properties of Lyapunov functions of multi-stable nonlinear systems but did not give the general method of how to construct them. Bhattacharya *et al.* [9] proposed a simple method to visualize a landscape in lieu of a vector field. But it represents neither Lyapunov stability nor transition rates among different attractors [9]. Ao and co-workers [10–14] suggested a general method to construct the quasi-potentials using symmetric and anti-symmetric decomposition. The method of Ao can be shown to be mathematically equivalent to our approach under certain conditions but it is not convenient for computing the transition rate between attractors in general. Xing [15] generalized Ao's work to prove that this approach is a mapping of a non-equilibrium dynamical system onto a Hamiltonian one, but misrepresented Ao's method as a generalized Helmholtz decomposition (see §3). Using stochastic systems where the dynamics manifests as the temporal evolution of the probability density function and is described by a Fokker–Planck equation, Wang *et al.* [16] championed the idea of a global potential function for gene networks using a decomposition that is similar but, again, not identical to the Helmholtz decomposition (see §3). Berglund [17] reviews Kramers' law and the computation of the transition rate in multi-stable dynamic system. He gives a rigorous mathematical discussion of the limits of applying Kramers' law in a non-gradient system [17].

Here we discuss several *ad hoc* methods used to construct or visualize a quasi-potential landscape and analyse whether they satisfy our first two conditions: (i) global metastability and (ii) agreement with the Freidlin–Wentzell theory for calculating the transition rates. Then we present the ‘normal decomposition’ as a systematic method that yields a quasi-potential that meets those conditions. Mathematical derivation and several examples are given to demonstrate the validity and the utility of our new approach.

## 3. Comparison of vector field decomposition

We first review several commonly used methods of vector field decompositions and the physical meaning of the components.

### 3.1. The Helmholtz decomposition

The *Hodge decomposition theorem* states that any sufficiently smooth (twice continuously differentiable), rapidly decaying vector field in *n*-dimensional space can be uniquely decomposed into the sum of a curl-free potential field d*x*/d*t* = **F**(**x**)*U*^{H}(*x*) and a divergence-free curl term [18]. Helmholtz's theorem is a special case of Hodge decomposition theorem in three dimensions, which can be written in the following equations:
3.1

The potential *U*^{H}(** x**) can be found by solving the following Poisson equations:
3.2

Because the curl term is divergence-free, i.e. (∇ × ** A**(

**)) =**

*x***0**, we have a Poisson equation, 3.3Without boundary conditions, the decomposition is not unique because any specific solution that satisfies equation (3.3) plus any harmonic function will also be qualified. But the decomposition is unique with the boundary conditions that

**F**(

**x**) →

**0**when

**approaches infinity. Therefore, the potential**

*x**U*

^{H}(

**) can be directly ‘read off’ from the dynamic equations (1.1) without a time-stepping process by solving the static Poisson equations with known right-hand terms and boundary conditions.**

*x*The meaning of this decomposition is that any dynamical field can be divided into two parts: a conservative irrotational field whose gradient connects the attractors or repellors (analogous to sinks and sources in the field theory) and a solenoidal curl field that is not influenced by sinks or sources. However, given that the divergence-free curl field has been used to explain the driving force of the gradient-free limit cycle [19], it is noteworthy that the curl field does not necessarily consist of closed trajectories. It can be made up of open trajectories with varying direction, as shown in figure 2. Nevertheless, this decomposition permits the study of the factors that influence the period of limit cycles (the curl term) and their stability (the gradient term) separately.

We can now determine whether the potential function obtained from the Helmholtz decomposition, *U*^{H}, satisfies the Lyapunov criterion for stability, our requirement (i), i.e. in particular, if d*U*^{H}/d*t* < 0 for ** x** ≠

***. For a two-dimensional system taking the time derivative of**

*x**U*

^{H}from equation (3.1), we obtain 3.4

Because we only know that ∇ · **F _{c}**(

**) = 0, (∇**

*x**U*

^{H},

**F**

_{c}) can be either positive or negative. Thus, d

*U*

^{H}/d

*t*is not necessarily non-positive, and hence Lyapunov stability is not guaranteed in this decomposition. Therefore,

*U*

^{H}does not satisfy our criterion (i) and it is not a measure of global metastability.

### 3.2. Decomposition based on the flux of the probability, *U* ∼ –ln*P*

An often taken, *ad hoc* approach to obtain a quasi-potential function *Ũ* and to visualize it as a landscape [16,20,21] stems from the intuitive notion that, in a probabilistic form of the system, the more stable states (lower potential *Ũ*) are more probable in the state space. Thus, the steady-state probability for each state *x* is then linked to a potential *U*^{prob}, according to the Boltzmann law
3.5

Again, because in general **F**(** x**) ≠ ∇

*U*

^{prob}, the driving force

**F**(

**) requires an additional component besides the gradient term. With the constraints of equation (3.5), i.e. a gradient force emanating from a probability flux as one component, Wang**

*x**et al*. [22] derived the following decomposition.

The time evolution of the probability density function *P*(** x**,

*t*) is governed by the Fokker–Planck equation, which can accurately describe the transition dynamics in a multi-stable dynamic system if the noise is Gaussian distribution, 3.6in which the drift term is the driving force

**F**(

*x*) of our system. Note that the diffusion term

*D*needs to be large enough for the probability density function

*P*(

**,**

*x**t*) to achieve a good coverage of the whole state space. The temporal change of the probability density can be written in terms of the probability flux

**J**(

**,**

*x**t*), 3.7

Because at steady states ∂*P*(** x**,

*t*)/∂

*t*= 0 the divergence of the steady-state probability flux

*J*_{ss}vanishes, i.e. ∇ ·

**J**(

_{ss}**,**

*x**t*) = 0, we obtain 3.8

In this vector field decomposition based on steady-state probability density, the driving force **F**(** x**) has been decomposed into two terms,
3.9Thus, the first term consists indeed of a gradient of a potential

*U*

^{prop}. The second term is the flux

*J*

_{ss}/

*P*

_{ss}and has a curl nature. However, it is not a curl in the same sense as in the Helmholtz decomposition because generally it is not divergence-free: ∇ · (

*J*

_{ss}/

*P*

_{ss})≠0. Only

*J*_{ss}itself is divergence-free. Thus this decomposition is not a Helmholtz decomposition as claimed elsewhere [21].

To test whether *U*^{prob} of this decomposition meets both our requirement (i) and our requirement (ii), we insert and **F** =− ∇*U*^{prob} + **F**_{c} into the Fokker–Planck equation,
3.10and obtain
3.11That is, the two vector fields are not perpendicular unless **F**_{c} = *J*_{ss}/*P*_{ss} is divergence-free—which is not necessarily so, as argued above. In other words, this decomposition, while motivated by *U*^{prob} ∼ −ln(*P*), is actually a field decomposition **F** =−∇*U* + **F**_{c} with the constraint (−∇*U*^{prob}, **F**_{c}) = *D*∇ · **F**_{c}.

We can now again determine whether the potential function obtained, *U*^{prob} ∼ −ln(*P*_{ss}), satisfies the Lyapunov criterion for stability, our requirement (i), i.e. in particular, if d*U*^{prob}/d*t* < 0 for ** x** ≠

***. For a two-dimensional system taking the time derivative of**

*x**U*

^{prob}, we obtain 3.12Thus, d

*U*

^{prob}/d

*t*is not necessarily non-positive, and hence Lyapunov stability is not guaranteed in this decomposition. Therefore, generally speaking,

*U*

^{prob}∼ −ln(

*P*

_{ss}) cannot be considered a measure of global metastability. One argument to save

*U*

^{prob}∼ −ln(

*P*) is that equation (3.11) will hold if the noise term

*D*goes to zero in the limit, i.e. if

*D*→ 0 then

*D*∇ ·

**F**

_{c}→ 0. It implies that the gradient component and the remainder

**F**

_{c}would be perpendicular to each other and that the Lyapunov criterion is also satisfied. However, the noise term

*D*cannot be too small. If

*D*indeed goes to zero, state transitions between attractors will never happen and

*P*

_{ss}will not be a global probability distribution but only a probability distribution in local attractors as a result of the balance between noise and driving forces.

We next examine the meaning for the case when the constraint is such that the gradient force and the remainder force are perpendicular to each other because this would guarantee satisfaction of the Lyapunov stability criterion.

### 3.3. The normal decomposition

From the discussion above, the constraint that the gradient term is perpendicular to the ‘remainder’ force **F**_{⊥} is of interest,
3.13

Taking the same approach as in equation (3.12), for a two-dimensional system, we obtain the time derivative for *U*^{norm} along the trajectory. Supposing that we have a Lyapunov function *U* for a two-dimensional system, its time derivative along any trajectory driven by the dynamic system in equation (1.1) is
3.14If the gradient of *U*^{norm} is normal to the remainder force **F**_{⊥}, i.e.
3.15then the function *U*^{norm} will decrease monotonically with time during the process before reaching the steady states, thus satisfying *Lyapunov's* condition for global metastability,
3.16We can therefore see that the condition for the normal decomposition has an obvious physical meaning, in that, for (∇*U*^{norm}, **F**_{⊥}) = 0, *U*^{norm} corresponds to a Lyapunov function *U*^{norm} of dynamical systems and can represent the global metastability. Thus, we seek a decomposition of a sufficiently smooth vector field into a conservative potential field *U*^{norm} and the remaining forces **F**_{⊥},
3.17

The physical interpretation is that for a ball in an attractor state that is perturbed to exit with least ‘energy’ against the system's ‘driving force’ **F**(** x**) that keeps it in the attractor, we can decompose the field such that

**F**

_{⊥}will not contribute to this process but only a gradient field

*U*

^{norm}contributes. On the basis of Freidlin–Wentzell's large deviation theory of a stochastic process discussed below [5],

*U*

^{norm}can thus be used to compute the transition rate in a non-equilibrium dynamic system.

The normal potential *U*^{norm} can also be read off the systems equations without a time-stepping solution. We can calculate the potential field *U*^{norm} as follows:
3.18

This equation is the Hamilton–Jacobi equation (HJE), which can be written in a component format
3.19The HJE provides for any non-equilibrium dynamic systems a general mathematical framework to construct a quasi-potential function that allows us to compare the relative stability of different attractors. This extends the global stability analysis of Lyapunov's theory to multi-stable dynamic systems. The quasi-potential function *U*^{norm} also epitomizes the epigenetic landscape that could explain the driving force behind cell differentiation in multicellular organisms, which can be used to calculate the transition rates between different cellular attractors (more details in §4 and §6)

Note that the existence of a solution for equation (3.19) is by no means guaranteed for all dynamic systems. However, here we deal only with a subset of problems encountered in biology rather than a general dynamic system. We now propose that the dynamical systems that represent organismal development admit such a decomposition. Further research is needed to find the conditions for the existence of such a decomposition and verify them in the real biological systems.

It is far from trivial to solve a nonlinear PDE such as the HJE, which is still a topic of active research in applied mathematics. Since there usually exist no analytical solutions, numerical methods with careful convergence and stability analysis have to be applied to solve the HJE in the real biological systems. A numerical method, the Newton–Raphson method, is used in this paper, which is described in the electronic supplementary material.

### 3.4 The relationship between the normal decomposition and the symmetric–antisymmetric decomposition

The symmetric–antisymmetric decomposition proposed by Ao *et al.* [23] is under certain conditions mathematically equivalent to the normal decomposition approach presented above.

#### 3.4.1. The symmetric–antisymmetric decomposition when **D** = *I*

**D**

For a dynamic system with governing equations in equation (1.1), the symmetric–antisymmetric decomposition can be achieved as follows [14]:
3.20where *U*^{ST} is the quasi-potential derived from the symmetric–antisymmetric decomposition*. S* is a symmetric operator,

**is an antisymmetric operator. Because it is usually difficult to find the analytical solutions for**

*T***and**

*S***,**

*T*

*S**can be solved as a*symmetric matrix and

**as an antisymmetric matrix in finite dimensional approximation.**

*T***is a diffusion matrix. Equation (1.1) can be rewritten as the following: 3.21Then**

*D*

**F**_{r}can be written as 3.22

*is an identity matrix.*

**I**To follow our perpendicular decomposition, the inner product of two terms is
3.23If ** D** is the identify matrix, i.e.

**= I, then, from equation (3.20), we have 3.24**

*D*And from equations (3.23) and (3.24), it is easy to show that (∇*U ^{ST}*,

*F*_{r}) = 0, i.e. the gradient of the quasi-potential ∇

*U*is perpendicular to

^{ST}

*F*_{r}. Therefore, the symmetric–antisymmetric decomposition is mathematically identical to our normal decomposition if the diffusion matrix

**is set to be an identity matrix**

*D***.**

*I*#### 3.4.2. Symmetric–antisymmetric decomposition when *D* ≠ *I*

*D*

The ensuing question is what the quasi-potential *U ^{ST}* obtained from the symmetric–antisymmetric decomposition will be if

**≠**

*D***. We can again evaluate this using the two conditions that we proposed in §2. First, we find that**

*I**U*decreases monotonically with time during the approach to the steady states, thus satisfying Lyapunov's condition for global metastability: 3.25We can therefore see whether the solution for equation (3.20) exists;

^{ST}*U*is indeed a Lyapunov function of a dynamical system that can represent global metastability.

^{ST}However, if ** D** is not an identity matrix

**, i.e. (∇**

*I**U*,

^{ST}

*F*_{r}) ≠ 0, the remainder

**F**

_{r}is not perpendicular to ∇

*U*and thus will contribute to the state transition. Then, on the basis of the Freidlin–Wentzell large deviation theory of a stochastic process (discussed below) [5],

^{ST}*U*cannot be used to compute the transition rate in a non-equilibrium dynamic system.

^{ST}## 4. Freidlin–Wentzell theory of large deviation in multi-stable systems and its relationship with the normal decomposition

### 4.1. Background

For a dynamic system with deterministic forces **F**(** x**,

*t*), 4.1Now let us consider that the system is under a stochastic perturbation

*ɛ*, 4.2If

*ɛ*is sufficiently small, the perturbed system will converge to the original dynamical system, i.e. ||

**−**

*x***|| → 0. However, if the perturbation is a random process with small average amplitude but with occasional large excursions, the perturbed dynamic system will behave differently. Freidlin & Wentzell [5] proposed a large deviation theory of stochastic process as a theoretical framework for the analysis of dynamical systems with multiple attractors. Supposing that a dynamical system satisfies the Langevin dynamics, the governing equations are described by the following ODEs: 4.3Supposing that a ball is perturbed to go from state by a stochastic process, one defines an action function**

*x*0*V*

*to measure the ‘energy’ barrier to be overcome for this transition, 4.4Here, the action function*

_{AB}*V*is defined as a time integral of the square of the ‘remainder’ of the dynamic equations (deviation from deterministic trajectory) over the whole trajectory

_{AB}*X*(

*t*) from attractor

**x**to

_{A}**x**. If a ball is only driven by the ‘forces’ specified in the deterministic part of the ODEs (equation 4.2), it will correspond to a ‘free fall’ in the ‘gravity field’ and the action is zero. If a ball is perturbed against the ‘forces’, the remainder term will not be zero. We integrate them over time along the whole trajectory and obtain the total action when the ball switches from attractor

_{B}**x**to

_{A}**x**. On the basis of the

_{B}*variation principle*, there exists a unique minimum action

*V*

**, which is an objective measure of the difficulty for the state transition in non-equilibrium dynamic systems.**

_{AB}### 4.2. The relationship between the Freidlin–Wentzell potential and normal decomposition

Although the Freidlin–Wentzell potential *V* and the normal potential *U*^{norm} are defined in different ways, they are mathematically related. For a dynamic system, **F** = − ∇*U*^{norm} + **F**_{⊥} with ( − ∇*U*^{norm}, **F**_{⊥}) = 0. We can rewrite the Wentzell potential as
4.5During the process of exiting an attractor, the LAP of a ball follows the governing equation:
4.6Intuitively, the LAP out of attractors will be driven by a landscape that reverses the original one, i.e. the hills are reverted to valleys, and valleys are changed to be hills. However, the reversal only holds for the gradient component because the normal remainder **F**_{⊥} still retain the same direction as before [24]. As long as we are within the same attractor, the Freidlin–Wentzell potential is exactly twice as big as the potential from the normal decomposition,

When a ‘ball’ transitions from one attractor to another, the Freidlin–Wentzell potential only accounts for the uphill ‘energy’, which is two times that of *U*^{norm}. Once it goes over the ‘saddle’ point (*S*_{1} and *S*_{2} in figure 1), the Wentzell potential *V* is zero for the remaining free-fall path. The same applies for a ball that transitions through *many* attractors in between: the Wentzell potential is equal to the sum of all uphill potentials between two points. All downhill paths contribute nothing to the Freidlin–Wentzell potential.

### 4.3. Wentzell potential and spontaneous transition rate

For a system that is under a large deviation perturbation: d*x*/d*t* = **F**(**x**,*t*) + *ɛ*·** ξ**(

*t*),

*V*, as defined in equation (4.4), the probability of spontaneous transition

_{AB}**x**

*→*

_{A}**x**

*, is [5] 4.7*

_{B}Accordingly, for the transition between *x** _{B}* →

*x**the probability is 4.8If*

_{A}*V*>

_{AB}*V*under random perturbation, a ball will stay longer in state

_{BA}**x**

_{B}, i.e. state

**x**

_{B}is globally more stable than state

**x**

_{A}and, thus, the directionality of the spontaneous transition points from state

**x**

_{A}to state

**x**

_{B}.

### 4.4. Wentzell potential and least action path of spontaneous transition

Although the Wentzell potential *V _{AB}* can be mathematically expressed as the normal potential

*U*

^{norm}it offers additional information that cannot be directly derived from

*U*

^{norm}. In determining

*V*, not only a minimum value is computed, but also a LAP

_{AB}*X*between two points is obtained.

_{AB}However, even if there is a LAP, a state transition **x**_{A} → **x**_{B} does not necessarily have to follow exactly this path only. Because the perturbation is random, the effective state transition can happen along any path that topologically connects these two states. The majority of probable paths fall near the area of the LAP, i.e. for any path **X _{t}**(

**x**

_{A}→

**x**

_{B}), the LAP is

**X**and 4.9

_{AB}Here, ‖ ‖ is defined as a norm, ‖ X ‖= It means that the bigger the noise, the further away the actual transition paths can deviate from the LAP **X _{AB}**.

## 5. Examples

In this section, two examples are presented to show which decomposition can give the correct transition rates between metastable attractors in nonlinear, high-dimensional dynamical systems. Quasi-potentials *U*(** x**) are computed using different decompositions: the Helmholtz decomposition (§3.1), the

*U*∼ –ln

*P*decomposition (§3.2) and the normal decomposition (§3.3). The transition barriers Δ

*U*are calculated based on the quasi-potentials constructed from these decompositions. The results are compared with the Wentzell action functions that represent the correct transition rates among different attractors.

### 5.1. Example 1. Spiral dynamic system with a single attractor

A two-dimensional non-gradient dynamical system reveals the effects of the non-vanishing remainder component of the vector field. The governing equations are 5.1The perpendicular decomposition of the vector field can be easily found from equations (3.18) and (3.19), 5.2

The quasi-potential functions and the remainders from the perpendicular decomposition are 5.3It can be easily verified that 5.4The Helmholtz decomposition can be constructed in a straightforward manner, 5.5

We have 5.6As mentioned before, the Helmholtz decomposition is not unique. The superposition of the function above with any harmonic function still satisfies the Helmholtz decomposition. However, this example illustrates the general characteristics of Helmholtz decomposition: the remainder and the gradient components are not necessarily perpendicular; thus, according to §4, it is not the right decomposition to give the correct transition rate.

Because *U*^{prob} in the –ln*P* decomposition (§3.2) depends on solving the steady states of the Fokker–Planck equation, which usually does not have an analytical solution, here we calculate *U* ∼ –ln*P* using a numerical method (finite difference method); see details in the electronic supplementary material. Different noise levels in the Fokker–Planck equation lead to different distributions *P*(**x**) in steady states. *U*^{prob} = −ln*P* ×*D* were calculated with different noise levels (*D* = 0.015, 0.045, 0.09) to cover a wider range of results from the –ln*P* decomposition.

The quasi-potentials at nine points (relative to the origin, i.e. *U*(0,0) = 0) from four different methods were calculated based on the methods described above and shown in figure 3. As can be recognized, the quasi-potential *U*^{H} (shown by circles) from the Helmholtz decomposition significantly departs from the Wentzell potential, and thus is not good for calculating the transition probability. The quasi-potential from the −ln*P* decomposition *U*^{Prob} has two problems. First, if the noise level *D* is small (*D* = 0.015), *U*^{prob} (shown by the symbol ‘x’) does not ‘cover’ the peaks on both sides. Second, even if the noise level *D* is large enough to sample the whole landscape, its quasi-potentials *U*^{prob} (shown by pentagons and diamonds) are not symmetrical owing to the asymmetric probability flux in the phase space. Here only the quasi-potential derived from normal decomposition *U*^{norm} (shown by squares) agrees with half of the values of the Freidlin–Wentzell potential (shown by triangles). The numerical method for calculating the Freidlin–Wentzell potential is in the electronic supplementary material.

### 5.2. Example 2. Two-dimensional dynamical system with four attractors

Because our original motivation was to determine the relative metastability of attractors in a multi-stable dynamical system, an example with four attractors is chosen to illustrate the various quasi-potential functions that can be extracted from the decomposition methods and allow us to compare the transition rates and global ordering of attractors they predict. (The Helmholtz decomposition is excluded from the analysis for its obvious invalidity.) Theoretical predictions from the two other decompositions are compared with the Freidlin–Wentzell action functions that can represent the transition rate between attractors caused by stochastic perturbation with large deviations.

The governing equations of our system are 5.7Because the dynamical system here is constructed to admit an analytical solution, it has a symmetric structure and can be decomposed exactly into a gradient part and a perpendicular remainder. However, the biological dynamic systems, such as the formulation of gene regulatory networks as rate equations for changes of gene expression, in general have no analytical solutions. Numerical methods are required to find the solutions approximately. (More details are included in the electronic supplementary material.) 5.8

The quasi-potential function *U*^{norm} from the normal decomposition is shown in figure 4*a*. It has four attractors designated as A, B, C and D, with the quasi-potential function *U*^{norm} as follows:
5.9

We can now examine the ordering of the attractors with respect to their relative metastability. Using equation (5.9), we find: *U*^{norm}(*D*) > *U*^{norm}(*B*) > *U*^{norm}(*C*) > *U*^{norm}(*A*). In this case, we have a transitive set, i.e. if *P _{D}*

_{→B}>

*P*

_{B}_{→D}and

*P*

_{B}_{→C}>

*P*

_{C}_{→B}then

*P*

_{D}_{→C}>

*P*

_{C}_{→D}and so forth. However, this does not need to be so in general (see §6, for more details).

We can verify that the quasi-potential function *U*^{norm}(*x*,*y*) in equation (5.9) satisfies the HJE (3.18). From the vector field plot in figure 4, it is clear that the gradient components contribute to the driving force toward the attractors (figure 4*b*) while the remaining perpendicular components capture the rotational movement in the phase space (figure 4*c*). However, because the rotational driving forces are perpendicular to the gradient components, they do not contribute to the transition between the different attractors.

Using the same method as in the previous example, the decomposition for extracting *U*^{prob} can be numerically computed with varying noise levels. In figure 4*d*, the results for the quasi-potential *U*^{prob} with various noise levels are plotted alongside *U*^{norm} from the normal decomposition. As shown in the previous example, the accuracy of *U*^{prob} suffers if the noise level is too small to cover the whole landscape. In this case, we need *D* = 20 to get a reasonably accurate landscape.

Then we need to calculate the Wentzell action function *V* as the gold standard to validate two decompositions here. Because the Wentzell action function is usually too complicated to be found analytically, equation (4.4) can be re-written in discrete form to find an approximate solution (see the electronic supplementary material for more details). We can calculate the LAP for the attractor transition from *x** _{A}* to

*x**and vice versa, and for the attractor transition from*

_{D}

*x**to*

_{A}

*x**and vice versa (figure 5). It is clear that the Wentzell action function applies during the ‘uphill’ process while it is zero during the ‘downhill’ process (figure 5). The LAP from attractor from*

_{C}

*x**to*

_{D}

*x**is different from the LAP from*

_{A}

*x**to*

_{A}

*x**. Here the Wentzell action function is smaller than for the transition from*

_{D}

*x**to*

_{A}

*x**because the attractor state*

_{D}

*x**has a ‘higher elevation’ than*

_{D}

*x**in the quasi-potential ‘landscape’*

_{A}*U*

^{norm}(figure 5

*a,b*). Note that the LAP

**x**→

_{A}**x**is also different from the reverse path

_{D}**x**→

_{D}**x**and the LAP

_{A}**x**→

_{A}**x**is different from the reverse path

_{C}**x**→

_{C}**x**(figure 5). Such ‘path irreversibility’ is a common characteristic of nonlinear dynamic systems.

_{A}On the basis of the Freidlin–Wentzell large deviation theory, the transition rate between two points can be determined from the transition barriers Δ*U* along the LAP. Instead of comparing the transition rates directly, Δ*U* calculated from the quasi-potential *U*^{norm} and by *U*^{prob} are compared with the Freidlin–Wentzell action functions *V*, as shown in figure 6. It is shown that *U*^{norm} agree well with the half of the Freidlin–Wentzell action functions *V* while *U*^{prob} mostly do not agree with it. It could be argued that the accuracy of *U*^{prob} could be further improved by using bigger noise and better numerical methods. However, from various numerical test, *U*_{prob} does not significantly improve its accuracy more than the example presented here.

## 6. Discussion

We show here that the quasi-potential *U*^{norm} obtained from the normal decomposition has a non-local meaning, in that the sign of Δ*U*^{norm} represents the relative stability of a pair of adjacent attractors with respect to the spontaneity to transition from one attractor to another. The spontaneity implies that the system can, in the presence of noise in the state variables ** x**, i.e. without an explicit external driving force that directs this transition as epitomized by signalling molecules in cell biology, move from a given state to the most probable neighbouring state. No matter what the transition path is between two neighbouring attractors,

**x**

_{A}and

**x**

_{B}, if it is always easier to transition from point

**x**to point

_{A}**x**and vice versa. By contrast, for the alternative decomposition based on

_{B}*U*

^{prob}∼ −ln

*P*, this equivalency between potential and transition rate is not exact although, for practical purposes, (where

*U*

_{s}is the potential at the ‘saddle’ (exit) point between these attractors) often scales with [16].

For the transition **x**_{A} → **x**_{C} in cases where the attractors **x**_{A} and **x**_{C} are not adjacent, but separated by attractor B that needs to be traversed, i.e. **x**_{A} → **x**_{B} → **x**_{C} with ‘saddle’ points between **x _{A}**,

**x**and saddle points between

_{B}**x**

_{A},

**x**, the above approach can be expanded [25], 6.1The case with more than one attractor in between can also be calculated with a similar approach to that in equation (6.1).

_{C}Theories of quasi-potential functions have been developed for the one-to-one relationship of two points (Freidlin–Wentzell theory, Kramers' law). Two caveats are in order here regarding the interpretation of Δ*U*^{norm} in systems with more than two attractors.

### 6.1. Δ*U* and ‘fate choice’ in multi-potent cells

First, the attractors with lower quasi-potentials *U*^{norm} have higher probability to be occupied when a system reaches its steady state—which may take a long time for a rugged landscape with small noise. This shall remind us that Δ*U*^{norm} alone does not allow one to predict the ‘fate choice’ of a stem cell in an attractor within a multi-attractor system that faces alternative transitions (representing its potential to differentiate into multiple distinct cell types spontaneously). The local landscape topography needs to be considered. Similar to the case for equilibrium systems, speaking of potentials, irrespective of whether path-dependent or not, one needs to bear in mind that Δ*U*^{norm} between two points is not a general measure of spontaneity of a process but that the latter depends on the local kinetic barriers, as illustrated in figure 1. Referring to the constellation in figure 1 it is clear that noise (or temperature) may be just high enough but not too high, such that one will find that the system in attractor A will prefer to move to B rather than to C (in the absence of catalysis). Only under sufficient noise (or high temperature) and time will the preference for C in steady state become manifest.

### 6.2. A global potential landscape?

Second, we now return to the more profound motivation behind the interest for a global quasi-potential function: the assessment of the ‘relative (meta)stability’ of attractor states, that is, their ordering with respect to spontaneous transitions into each other. Such a global behaviour of a system, as opposed to local, peri-attractor behaviour, is the scale of system change at which complex biological processes, such as cell differentiation during the development, take place. But is such an ordering possible? Can we directly compare the values of the quasi-potentials ) for *m* attractors and extract meaningful information as to long-term fate choice—aside from the aforementioned issue of kinetic inhibition? This is tempting to ask because, unlike the Freidlin–Wentzell potential difference *ΔV*, which is defined for pairs of states, we can compute from the system equations directly the potential value *U*^{norm} (*S*) for any state *S*. If *U*^{norm} captures the global behaviour, one would postulate that our conjectured third condition, the transitivity of any three points *a, b, c* with respect to their quasi-potential *U*, would be satisfied: if *P _{a}*

_{→b}>

*P*

_{b}_{→a}and

*P*

_{b}_{→c}>

*P*

_{c}_{→b}, then

*P*

_{a}_{→c}>

*P*

_{c}_{→a}. Example 2 constitutes such a transitive system. However, there is no guarantee that transitivity is satisfied for all systems. Although Δ

*U*

^{norm}scales monotonically with Δ

*V*(Wentzell action function), when it is defined for points within the same basin of attraction (most usefully, as the difference between the attractor state and the exit state at the saddle), there is no simple mathematical relationship between Δ

*U*

^{norm}and Δ

*V*along the transition trajectory between two attractors and it thus is ‘agnostic’ for global ordering.

Only in those cases where there is transitivity among all the potential values of attractors *i*, (consistency) in the directed graph in which the edge arrows depict the direction of higher transition probability *P*(**x*** _{i}*→

**x**

*) between pairs of attractors*

_{j}**x**

*and*

_{i}**x**

*, can we speak of a ‘global’ landscape. This condition surely will not be satisfied for systems exhibiting biological phenomena such as circadian cycle, cell cycle and other biological limit cycles. However, in some classes of systems, as lucidly epitomized in the qualitative picture of Waddington's epigenetic landscape that describes the essential features of cell development or in the fitness landscape of evolution, such global landscapes appear to adequately capture the developmental or evolutionary dynamics and constraints. In differentiation, it visualizes the barriers of trans-differentiation or reprogramming from one cell type to another. This suggests the possibility of a global ordering of (meta)stability for a particular yet to be defined subclass of networks to which our example 2 belongs, in which the global transitivity of*

_{j}*U*is satisfied. Thus, it could be that evolved networks that govern acyclic phenotype transitions, such as development or evolution, have for some unknown reasons a quasi-potential landscape that satisfies the transitivity condition (iii)—which is a subset of multi-attractor systems. The facility to compute for each attractor

*i*of a biological network a quasi-potential

*U*opens the possibility to study this phenomenon.

_{i}### 6.3. Present and future biological significance

Knowledge of the quasi-potential *U*^{norm}, if it can be computed, has direct significance for understanding and manipulating regulation of developmental paths in multicellular systems. This would require an accurate dynamic system description of the relevant gene regulatory network and measurement of the state vector for an attractor state of interest—corresponding to biologically defined cell types. can then be computed from the normal decomposition of the system equations for the attractor states , which currently can be approximately measured as transcriptome or proteome. in turn will afford a ‘global ordering’ of the cell types *i* with regard to their metastability. The ordering of a set of different cell types can be used to predict their relative stability in terms of the ‘ease’ (probability) of transitions into ‘neighbouring’ cell types, i.e. transitions between attractor states under noise or perturbation. The gradient of *U*^{norm} can be interpreted as the driving force that imposes directionality to the process of cell differentiation. This ultimately can help explain how the gene–gene interactions encoded in the genomic DNA sequence unfold to manifest as a robust, spontaneous and directed developmental process.

More concretely, because *U*^{norm} of two different cell types also provide an objective measure to compare their relative (meta)stability, we could in principle find the best subset of genes whose expression needs to be jointly perturbed to minimize the ‘transition barrier’ between these two cellular attractor states. This could pave the way towards a rational design of cell reprogramming protocols for which it is currently only known that combinations of a number of genes need to be perturbed, the determination of which however is purely empirical.

At present, the major limitation to realizing the above utility of the quasi-potential formalism is that information about the structure of the gene regulatory network, ** F**(x), is only sketchy at best, limiting its description as a dynamic system (equation (1.1))—the prerequisite for determining

*U*

^{norm}(x). As we learn more details about

*(*

**F***x*) the quasi-potential will gain in importance accordingly.

## Acknowledgements

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), Canadian Institutes of Health Research (CIHR) and Alberta Innovates, as well as by the Institute for Systems Biology (ISB), Seattle.

- Received May 28, 2012.
- Accepted August 3, 2012.

- This journal is © 2012 The Royal Society