Analysis of stochastic bifurcations with phase portraits

We propose a method to obtain phase portraits for stochastic systems. Starting from the Fokker-Planck equation, we separate the dynamics into a convective and a diffusive part. We show that stable and unstable fixed points of the convective field correspond to maxima and minima of the stationary probability distribution if the probability current vanishes at these points. Stochastic phase portraits, which are vector plots of the convective field, therefore indicate the extrema of the stationary distribution and can be used to identify stochastic bifurcations that change the number and stability of these extrema. We show that limit cycles in stochastic phase portraits can indicate ridges of the probability distribution, and we identify a novel type of stochastic bifurcation, where the probability maximum moves to the edge of the system through a gap between the two nullclines of the convective field.


Introduction
Networks of interacting species play an important role in many disciplines like for example chemistry, ecology, or systems biology [1]. A powerful modeling approach for such networks are ODE (ordinary differential equation) models, which are often written in the form of a dynamical system, The functionf ðxÞ is called the deterministic drift and has the mathematical structure of a vector field [2]. In particular in one-and two-dimensional systems, the general dynamical behavior can be directly deduced from vector plots off ðxÞ [3]. Based on these, the qualitatively different types of trajectories and the invariant manifolds (fixed points, limit cycles) can be graphically represented in the form of so-called phase portraits [4]. When control parameters, such as reaction rates, are changed, dynamical systems can undergo bifurcations, during which fixed points or limit cycles are destroyed or generated, or they change their stability. For instance, when the rate of predation in a predator-prey system is decreased, eventually a point will be reached where the predator can no longer find enough food and becomes extinct. Local bifurcations can be calculated fromf ðxÞ by determining the PLOS  fixed points and their stability in dependence of the parameters. Global bifurcations, which destroy or generate limit cycles of a finite size, are much more difficult to deduce fromf ðxÞ. In order to obtain a general first overview over the different types of bifurcations that occur in the system, it is useful to plot phase portraits for different values of the control parameters and to compare them to each other.
The above deterministic and continuous description is an idealization that is often too simple [5][6][7]. ODE models are based on mean reaction rates and are a good approximation when concentrations are large and systems are well mixed. They have been applied successfully to many different biological and chemical networks [8][9][10]. However, species concentrations cannot be treated as continuous quantities when numbers become small, as there are no half molecules, proteins, or animals. In this case, one should use a discrete state space. Furthermore, the randomness of reactions can lead to considerable deviations from the deterministic dynamics [11][12][13]. For instance, it can lead to accidental extinctions of species, to the appearance of quasi-cycles in predator-prey models [14], or to switching between two types of behavior, as in the foraging efforts of ant colonies [15]. Additionally, the same deterministic equation can arise from different microscopic models, either by using a reaction that produces one molecule often, or a reaction that produces several molecules at once, but less often. Such a bursty reaction occurs for example during translation in protein synthesis [13,16,17]. Different burst sizes lead to a different intrinsic noise of the system: Not only is the variance of the noise distribution proportional to the burst size, but the noise becomes asymmetric when only production but not the destruction of molecules is bursty. Such an asymmetry can shift the mean of the stationary distribution with respect to the value expected from the ODE model. [5,18] In order to assess the relevance of stochasticity, it is useful to express the variables in terms of the system size or reactor volume N, which is a measure for the total number of reactants in the system. It is well known that the fluctuations in a stochastic system scale only with ffiffiffiffi N p , whereas the mean values scale with N. This means that the signal-to-noise ratio decreases with N. In the limit N ! 1, the ODE description is recovered. [5] There exist several modeling techniques for such stochastic systems, which describe the system on different levels of resolution. The most important of these techniques is the Master equation, which models the time evolution of the numbers of the different molecules as a Markov process with a discrete state space and continuous time [6]. As the time evolution depends only on the present state and not on the past, memory effects are ignored. Furthermore, a wellmixed system is assumed since no spatial effects are considered. Despite its simplicity, the Master equation cannot usually be solved analytically. For that reason, various approximations are used, which simplify the calculations but at the same time preserve important features of the stochastic description.
One of these approximations is the Fokker-Planck equation, which is obtained from the Master equation by keeping the lowest-order terms of the so-called Kramers-Moyal expansion [19]. This expansion uses a coarse-grained, continuous state space and is essentially a Taylor expansion to second order in the changes of these continuous state variables. The Fokker-Planck equation becomes more accurate when more molecules are in the system. It works however still astonishingly well for systems with only a few tenths of molecules, even when ODE descriptions are already very bad [20]. The Fokker-Planck equation will also be the starting point for our investigation. A competing approach is the system-size expansion [5], which uses the inverse system size as an expansion parameter. It is often considered to be more systematic than the Fokker-Planck equation, but nevertheless frequently leads to results which are worse, at least when carried out only to the first two terms [20].
Besides those generic approaches, a vast amount of more specific techniques have been developed, which allow for all kinds of analyses of stochastic systems, such as effective deterministic descriptions [21], most probable dynamics [22,23], stability analysis [24], mean switching times of multistable systems [6], stochastic resonant cycles [14], or piecewise deterministic Markov processes [17,25].
So far, these techniques were rarely used to investigate the effect of stochasticity on the bifurcations of reaction networks. When the number of molecules is large and the dynamics resemble that of the deterministic ODE system, attractors of the dynamical system give rise to sharp maxima of the stationary probability density of the stochastic system. For this reason, bifurcations of stochastic systems are often understood as qualitative changes in the topology of the maxima and ridges of the stationary distribution. Those maxima can merge or split, such that a multistable distribution becomes monostable or vice versa. This class of bifurcations is called phenomenological bifurcations, or p-bifurcations [26]. We have investigated a stochastic bifurcation of this type in an analysis [27] of the Schloegl model [28] where we showed that burst noise shifts the parameter values of the saddle-node bifurcation. Burst noise can thus lead to bistability in a parameter region where the deterministic model or the model with conventional noise is monostable, and vice versa. Similar effects were found in other types of systems, such as in noisy bistable fractional-order systems [29], or for the Duffing-Van der Pol oscillator [30,31].
Besides p-bifurcations, Arnold [26] defines dynamical bifurcations (d-bifurcations). Dynamical bifurcations correspond to a separation of dynamics into state space regions between which no transitions are possible. This definition is chosen in analogy to bifurcations in ODE systems. For this reason, several subsequent publications [32,33] concentrate mainly on d-bifurcations. However, these bifurcations can occur only under the condition that the intrinsic noise of the system vanishes at the boundary between the different dynamical regions. This cannot happen for the kind of chemical reaction networks that we are interested in.
Therefore, we focus on p-bifurcations, which are very common in stochastic reaction networks. Since the stationary distribution of the probability density contains no information about the dynamics of the system, our study will not be confined to the calculation of stationary densities. Instead, we want to suggest a stochastic analogue of phase portraits, which convey a better understanding of the dynamics of the system and give an immediate insight into the possible bifurcations. Such a stochastic phase portrait has been used earlier by H. Qian at al. [34] in an investigation of a chemical oscillator with intrinsic noise in order to demonstrate the transition from a steady state to an oscillating state. Our work formalizes and generalizes this approach.
A related idea is pursued in [22], where the "most probable" trajectories are calculated, starting from an initial delta distribution. In contrast, our method gives information on the probability flows everywhere in state space starting from a constant distribution and is therefore not limited to unimodal distributions. This is, of course, crucial to the investigation of stochastic bifurcations. We will discuss the connection between these two formalisms further in the final section.
This article is structured as follows: Section 1 summarizes the methods that will be used, introduces stochastic phase portraits and shows that extrema of the stationary probability distribution coincide with fixed points of the convective field when the probability current vanishes.
The second part of the paper gives several examples for the application of these phase portraits, demonstrating the power of the method. Thus, the bifurcation diagram of a model for foraging ant colonies [15] will be analyzed without need to solve the Fokker-Planck equation, and the relation between limit cycles of the convective field and ridges of the probability distribution will be explored using the Rosenzweig-MacArthur predator-prey model.
As part of this investigation, we will encounter a special type of bifurcation that can only emerge in stochastic systems and has no deterministic counterpart. Finally, we will summarize and discuss our results in section 3.

Chemical reaction networks
We consider reaction networks for species X i that undergo a set of reactions with so-called stoichiometric constants σ ij and ρ ij and reaction rates μ j .
To compactify the notation, one defines [5] the stoichiometric matrix and the propensity vector which equals the reaction rate μ j per molecule(s) times the probability of the involved molecules to meet. In this notation n i is the molecule number of species X i and N is the system size or reactor volume, which can be interpreted as an inverse measure of the intrinsic noise of the system. With these definitions, the set of reactions (2) can be written in form of the so-called (chemical) master equation, with the step operator E i , which is defined via for an arbitrary function f and an integer α.

The multidimensional Fokker-Planck equation
The master equation yields a time-continuous, discrete-state Markov model of the system. As mentioned above, the analysis of the dynamics becomes simpler when using a continuous state space. Therefore we introduce the species concentrations x i = n i /N and base our analysis on the multidimensional Fokker-Planck equation (FPE) [6] @pðx; tÞ @t ¼ À Here, we have introduced the deterministic drift and the diffusion matrix From the definition of (positive) definiteness hx, A Á xi > 0, it follows that D is positive definite.
The drift term alone determines the ODE description of the system, which becomes exact in the limit of infinitely many molecules, The FPE has the form of a continuity equation, @p @t ¼ Àr Áj with the probability current Definingã ðxÞ ¼f ðxÞ À 1 2N with i ¼x i =jx i j being the unity vector in direction x i , the current can be written as The first termã i ðxÞpðx; tÞ is usually called the convection term, the second one the diffusion term [35]. To avoid confusion with the driftf and diffusion D terms in the FPE, we will refer to these two contributions to the probability current as the convective and diffusive current.
The expression for the convective currentj c ¼ pã has the same form as that for the electric current,j ¼ sẼ, and we therefore callã the convective field. This field depends only onf and D, but not on the probability density pðx; tÞ. Its effect is to pump probability towards the preferred states of the system, just like an electric field pushes positive charges towards its sinks.
In absence of the diffusive current, the convective current would converge to the attractors of the vector fieldã, producing a stationary probability distribution that is a (properly scaled) delta distribution on the attractors. The diffusive current, on the other hand, is directed toward decreasing p and thus tends to flatten steep slopes, ironing out sharp initial peaks and giving rise to a finite width of the stationary distribution. When a constant probability density is chosen as initial distribution, the diffusive current vanishes, and the convective field gives the direction of the current during the first infinitesimal time interval.
The stationary distribution is given by the condition @pðx;tÞ @t ¼ 0, which meansr Áj ¼ 0. In a one-dimensional closed system, wherej must vanish at the boundaries, this means thatj ¼ 0 in the stationary state. The convective and diffusive currents cancel each other. This means that the stationary state of a one-dimensional closed system satisfies detailed balance. Fig 1  illustrates the effect of the two fields on a narrow initial distribution for a one-dimensional system. The initial peak moves towards the stationary solution and becomes broader until in the stationary state the two currents cancel each other. In more than one dimension, stationary solutions without detailed balance (j 6 ¼ 0) are also possible, as the conditionr Áj ¼ 0 can be satisfied with nonzero currents. Such states are often called non-equilibrium steady states. The conditionr Áj ¼ 0 implies for these states thatj ¼r ÂÃ with a vector fieldÃðxÞ, i.e.,j is solenoidal.

Favorable and unfavorable states
We define favorable states as maxima of the stationary probability density at which the probability currentj vanishes. At these maxima, we have @p @x k ¼ 0 8k in (13). The conditionj ¼ 0 applied to (13) then givesã ¼ 0. This means that favorable states are always fixed points ofã. However the opposite needs not to be true, which we will discuss later.
This definition of favorable states is a generalization of the condition which gives the maxima and minima of the FPE in one-dimensional systems [7]. The calculation of the favorable states can be done without explicitly solving the FPE. In fact, this calculation is mathematically almost identical to the calculation of the fixed points in the deterministic description, which are obtained fromf ðxÞ ¼ 0. The only difference is thatf is replaced withã in the stochastic description.
In analogy to unstable fixed points of deterministic systems, we also consider minima of the stationary distribution p 1 ðxÞ at which the probability current vanishes. These are also fixed points ofã, and we call them unfavorable states of the system.
Favorable and unfavorable states are associated with a different stability of the fixed point ofã: Setting (13) equal to zero and differentiating with respect to a component x m , we obtain Using @p This can be rewritten as a matrix equation Since D is positive definite, its inverse is also positive definite. The multiplication with a positive definite matrix does not change the definiteness of a matrix. This means that the Hesse matrix h is positive definite when J s is positive definite, and negative definite when J s is negative definite. It follows that favorable (unfavorable) states are stable (unstable) fixed points ofã. The classification of these states can therefore be done with the same type of calculation as in the deterministic model, only withf being replaced withã.
There are also maximax Ã of the stationary probability density at whichrp 1 ðx Ã Þ ¼ 0 but jðx Ã Þ 6 ¼ 0. These points are not fixed points ofã. An example is given in Fig 2. If the deterministic system has a stable limit cycle with an unstable spiral in the center, the stochastic system can be expected to have an unfavorable state sitting in a crater surrounded by a ridge. Typically, the current along this ridge will not be constant, being faster at saddle points and slower at local maxima of the stationary probability density. These maxima cannot be classified as favorable states of the system, as trajectories are not attracted to them but pass through them. They are rather slow transient states.
The reverse situation, thatãðx 0 Þ ¼ 0 but thatx 0 is no extremum of p 1 can also be imagined. The current density at such a point must be different from zero. Although we never encountered such a "false positive" in the models we investigated, we could so far not show under which conditions they may occur or whether they can be ruled out completely. We will therefore always compare the calculated fixed points ofã with simulation results of the stochastic system, to verify that such a fixed point indeed corresponds to a favorable state of the system.
The authors of [34] discuss this question also, and they argue (in Appendix III) that the three conditionsrp ¼ 0,ã ¼ 0, andj ¼ 0 coincide. While this might be correct in some special cases, they do not give necessary and sufficient conditions that are generally valid. As mentioned above, slow transient states are a counterexample whererp ¼ 0, butã andj are different from 0.

Stochastic phase portraits
We will explore further below the usefulness of what we want to call stochastic phase portraits (SPPs), which show the qualitatively different types of trajectories defined by the vector field aðxÞ. This means that we consider the dynamical system for which we draw (conventional) phase portraits. As shown above, the stable fixed points of (18) coincide with the favorable states of the underlying stochastic system when the probability current vanishes at these points. Furthermore, we will use the trajectories calculated by solving (18) numerically as an indicator of the probability flow in the system. This will be especially useful to determine the occurrence and position of stochastic limit cycles. As these trajectories include only the convection flow of the probability current and are completely ignorant of the initial conditions of the stochastic system, they can, however, provide no information about the time evolution of a stochastic system with a specific initial distribution. In particular, one must not confuse these trajectories with sample paths of the stochastic system, nor with the evolution of its mean or of its maximum. The latter was investigated in a study of the most probable dynamics in [22].

Example: Foraging colony
In [15] the authors propose a model for the emergence of bistability in the foraging behavior of an ant colony, which we want to use as a first application of the SPPs. The model introduces two ant "species" X 1 and X 2 that feed on two different food sources. There are two types of "reactions". The first one corresponds to the recruiting of an ant of one species by the other species, so that it changes its feeding source. The other type of reaction is the spontaneous switching of food source by an ant. The authors in [15] find that the system exhibits either bistable behavior, where most ants visit the same food source and switch unregularly and within a short time to the other source or a monostable behavior where both food sources are frequented by approximately half the ants. Whether mono-or bistable behavior occurs depends on the quantity with N C ) N implying bistability and N C ( N monostability. In particular, monostable behavior is the only type of behavior found in the deterministic limit of very large colony sizes N. Since the total number of ants is constant, the number of ants of species 2 is given by n 2 = N − n 1 . This is equivalent to the condition for the species concentrations x 1 + x 2 = 1. As a result of this conservation law, the FPE of this system can be reduced to one-dimension. The authors solved it analytically in the stationary case, which leads to the probability distributions shown in Fig 3. Instead of solving the FPE explicitly, we will calculate the favorable states directly using (12) and show the SPP of the system. Even though the model is effectively one-dimensional, we will always use both variables, x 1 and x 2 .
Drift and diffusion are given bỹ and D ¼ This leads to a convection flow of The roots of (24) are x 1 = x 2 = 1/2. This means that we only expect one extremum of the stationary probability distribution, which lies at the same position as the deterministic fixed point that can be calculated from (22). However, this does not mean that the stochastic behavior of the system is similar to the deterministic one: The Jacobian of the deterministic system (22) is with eigenvalues λ 2 {0, −2}. The deterministic attractor of the system is therefore a stable fixed point line at x 1 = x 2 = 1/2.
The stochastic Jacobian, which is based onã, reads though with eigenvalues λ 2 {0, −2( − r/N)}. The stochastic fixed point at x 1 = x 2 = 1/2 is only stable for r/N < . Otherwise, the fixed point is a minimum, and the maxima of the probability distribution are located at the boundary. This is exactly the same result as obtained from the full solution of the FPE in (21). The fact that favorable states (or maxima of the probability distribution) can also emerge at the boundary of the state space, even though the deterministic drift at this points is not zero, can be interpreted as a special feature of stochastic systems which cannot occur in a deterministic description. We will discuss this in more detail in our second example. Of course, our findings so far provide no new insights into the ant model compared to the original investigation in [15]. We want to note though that our results can be obtained without the need to solve any differential equation and by following a formalism that is very similar to the deterministic analysis of the model.
Based on (24) and the stability of the fixed point from (26), we are now able to draw the SPP for the ant model, which is depicted in Fig 4.

General properties.
Next, we focus on a truly two-dimensional model that can show limit cycles in the deterministic version. This is the so-called Rosenzweig-MacArthur model for predator-prey interactions [9]. It is usually given in terms of ODEs, but it can be also represented by the following set of reactions: prey birth predator À prey interaction The immigration reactions are only introduced to prevent spontaneous extinctions (the socalled Keizer's paradox [36]) and can be neglected for the deterministic analysis. Drift vector and diffusion matrix for this system are easy to calculate and read Depending on the effective parameters κ = Ac/d and ¼ dc bd , the deterministic system corresponding to (28) can either exhibit stable populations of both predator and prey, a limit cycle where the populations of both species oscillate, or the extinction of the predator [37,38]. A phase diagram showing the parameter dependency of these regimes is depicted in Fig 5, along with the (deterministic) phase portraits for two sets of parameter values in Figs 6 and 7. Analysis of stochastic bifurcations with phase portraits 2.2.2 The effect of stochasticity when the deterministic system has a stable fixed point. Fig 8 shows SPPs for different system sizes N with reaction rates chosen such that the deterministic system has a stable fixed point. In order to assess the significance of the SPPs, we also performed Gillespie simulations of the corresponding reaction system using the free software tool Dizzy Gillespie [39]. The time series data from these simulations was used to generate density histograms, which correspond to the respective stationary probability distributions. These are shown as yellow clouds overlayed over the SPPs. For the large system size N = 1000, the SPP shown in Fig 8(a) agrees with the deterministic phase portrait in Fig 6. Fluctuations around the maximum value, which are obtained from the stochastic simulation, are rather small. For the smaller system size of N = 80, we obtain a similar SPP, however, the probability distribution obtained from the stochastic simulation is broadened due to the increasing noise. The maximum of the probability distribution coincides well with the stable fixed point of the SPP. For N = 50, the SPP remains unchanged, but the stochastic data become qualitatively different: While there is only one density maximum at approximately (3, 0.7) for N = 80, there emerges a second maximum at the boundary of the phase space at approximately (4, 0) for N = 50. We can interpret this emergence of a second maximum as a stochastic bifurcation, which we will discuss in more detail in the next subsection. For N = 35, the boundary maximum becomes more important, while the central maximum is almost gone. At N = 10, a second bifurcation has occurred, which is directly visible in the SPP: As the two nullclines cease to intersect, the central maximum vanishes and all probability becomes concentrated close to the boundary maximum at x = 4.
The SPP for N = 10 provides also a very clear understanding of how the maximum at the boundary emerges: As no fixed point is left inside the phase space the probability follows the convective field until a boundary (in this case the x-axis) is reached. From this position, the probability flow can't follow the convective field anymore since negative values must retain probability zero. Therefore, the probability flow slides along the x-axis until the x-component of the convective field equals zero. Naturally, this happens at the position where the x-nullcline intersects the x-axis. The boundary maxima emerge therefore always at the intersection points of a nullcline with its corresponding axis.

Nullcline gap bifurcation.
With this understanding of the boundary maxima, we reconsider the bifurcation that leads to their formation, which occurs for our system between N = 80 and N = 50. Comparing the SPPs for those two system sizes (Fig 8(b) and 8(c)), one sees that the two predator nullclines at the bottom exchange their direction. This is also illustrated in the exaggerated scheme in Fig 9. The dashed lines indicate the y-nullclines. As their    (N ( N C , Fig (b)) or is caught inside (N ) N C , Fig (a)) so that no boundary maxima can emerge. https://doi.org/10.1371/journal.pone.0196126.g009 Analysis of stochastic bifurcations with phase portraits behavior close to the x-axis changes qualitatively, the direction of the convective flow also changes: When the convection vectors in the gap point upwards (Fig 9a), the probability moves away from the x-axis and cannot produce a boundary maximum. When the vectors point downwards (Fig 9b), the probability flow in the gap region moves towards the axis, and a boundary maximum will emerge.
Since the formation of a boundary maximum is directly coupled to the change of the nullcline gap, we call this type of bifurcation a nullcline gap bifurcation.

2.2.4
The effect of stochasticity when the deterministic system has a limit cycle. Fig 10  shows SPPs for the parameter set that yields a limit cycle in the deterministic model.
Again, there is a good agreement between the SPP shown in Fig 10(a) and the deterministic phase portrait shown in Fig 7. The stationary probability obtained from the stochastic simulation is concentrated at the location of the limit cycle of the SPP. A closer look reveals though that some states have a much higher probability than others (indicated by a darker yellow coloring). These are the slow transient states discussed in section 3, where p 1 has maxima with non-vanishing probability current.
For N = 80, the probability distribution again becomes broader. Around N = 65 the Nullcline Gap bifurcation occurs, leading to a boundary maximum at x % 9 for the smaller system size of N = 50.

Discussion
We have presented a method to obtain phase portraits for stochastic systems (SPPs) as vector plots of the convective field obtained from the Fokker-Planck equation (FPE). This field can be obtained directly from the drift and diffusion terms without the need to solve the FPE. We showed that stable (unstable) fixed points of the convective field correspond to maxima (minima) of the stationary probability distribution if the probability flow vanishes at this point. This means that for a one-dimensional system the SPPs always reproduce correctly the maxima and minima of the stationary probability distribution (as obtained from the Fokker-Planck equation). We demonstrated this using a model of foraging ants [15], for which the SPP shows the change of stability of the symmetric fixed point as the system size is changed. Earlier, an explicit solution of the FPE was required in order to derive this result.
Furthermore, we considered a predator-prey model [9], for which the stationary probability flow does not vanish. We could also verify that the SPPs yield an accurate description of the stochastic dynamics of the system, which were obtained from Gillespie simulations. In particular, we found that boundary maxima of the stationary probability distribution can also be predicted from the SPPs. The emergence of these boundary maxima happens as a special type of bifurcation which can only appear in stochastic systems, and which we called nullcline-gap bifurcation.
These successes of SPPs suggest that they are for many systems a good indicator of where the probability concentrates in the stationary state. However, we were so far not able to determine whether a fixed point of the convective fieldã always leads to an extremum in the stationary probability distribution for multidimensional systems. So far, we did not yet find a counterexample. Similarly, it is unclear whether limit cycles in the SPPs correspond always to a non-equilibrium steady-state withj 6 ¼ 0 of the underlying stochastic system. Furthermore, although our example showed a higher concentration of the probability density along the limit cycle, we could not prove that limit cycles in the SPP are always associated with a crater-shaped stationary probability distribution, and if so, whether the ridge of this crater will always coincide with the course of a limit cycle in the SPP. In any case, as the SPPs are based on the Fokker-Planck equation and not on the Master equation, the agreement of the SPPs with the results of stochastic simulations cannot be better than that of the stationary distribution resulting from Fokker-Planck equation with the true stationary distribution. This indicates a general limitation of our method: While the Fokker-Planck equation is able to predict many features of stochastic systems -like for example noise-induced bistability [15,27]-very accurately when molecule numbers are not too small, the FPE breaks down for other systems. This happens in particular when stochastic effects are directly induced by the discreteness of the involved particle numbers [40]. Clearly it would be useless to apply the methods outlined in this paper to such systems. See also [20,41] for an overview of the limitations of the Fokker-Planck equation.
The results obtained so far suggest further directions of study. In particular, it would be helpful if one could extract information about the relative height of the maxima of the stationary probability distribution from the SPPs. We have seen examples where the SPP suggests a maximum at the boundary and one inside the state space, while in reality the boundary maximum was so high that it was almost impossible to notice the small maximum inside the state space. The size of the maxima might be correlated with the relative sizes of the basins of attraction or with the density of the probability flow of the convective field. Similarly, the SPPs give so far no indication of the widths of the respective maxima of the system. We suspect that progress in this direction can be obtained by exploiting the properties of the diffusion matrix.
Finally, we address the question how our method relates to the work of Cheng et. al. [22], where the most probable trajectories of a stochastic system are investigated. The authors start from a sharp initial probability distribution and calculate analytically the time evolution of the maximum. However, this calculation is based on the premise that the probability distribution shows exactly one maximum at all times t. In contrast, our approach has no such restrictions and describes the stationary maxima and minima also precisely. But our approach does not describe transient dynamics, except for the initial dynamics that occur when starting with a flat initial probability distribution. The two approaches are thus complementary to each other. Together, they provide an analytical means to investigate the extrema of multistable Fokker-Planck equations, which are otherwise only solvable by numerical methods.