Basin stability in delayed dynamics

Basin stability (BS) is a universal concept for complex systems studies, which focuses on the volume of the basin of attraction instead of the traditional linearization-based approach. It has a lot of applications in real-world systems especially in dynamical systems with a phenomenon of multi-stability, which is even more ubiquitous in delayed dynamics such as the firing neurons, the climatological processes, and the power grids. Due to the infinite dimensional property of the space for the initial values, how to properly define the basin’s volume for delayed dynamics remains a fundamental problem. We propose here a technique which projects the infinite dimensional initial state space to a finite-dimensional Euclidean space by expanding the initial function along with different orthogonal or nonorthogonal basis. A generalized concept of basin’s volume in delayed dynamics and a highly practicable calculating algorithm with a cross-validation procedure are provided to numerically estimate the basin of attraction in delayed dynamics. We show potential applicabilities of this approach by applying it to study several representative systems of biological or/and physical significance, including the delayed Hopfield neuronal model with multistability and delayed complex networks with synchronization dynamics.


I. FUNCTION SPACE FOR INITIAL VALUES
The Bernstein basis. In addition to the two bases introduced in the main text to construct the function space for the initial values, the Bernstein basis, which is the third basis we use in the main text, can be expressed as B j,m (t) = C j m t j (1 − t) m− j , j = 0, 1, · · · , m, 0 ≤ t ≤ 1, where C j m is the binomial coefficient [1]. The Bernstein basis, though nonorthogonal, offers valuable insight into its geometrical behavior and can be used to approximate any continuous function f on [0, 1], that is, where the limit is interpreted in a sense of uniform convergence for continuous functions [2]. For any given n, the Bernstein basis is changing, which is different from the Fourier series where a priori given is an orthogonal basis that contains an infinite number of basis functions. Moreover, the coefficients in the above approximation correspond to the values of f on the lattice, so that the Bernstein basis is always used in the curve or surface interpolation. As for our problem, for any given n, we set f i in C (n,α) as f i (x) = B i,n (x) and let all coefficients a i satisfy n i=1 |a i | < α for some given α. Then, akin to the algorithms proposed in the main text, we can compute the nth-order approximate basin stability (BS) for delayed dynamics.
As a matter of fact, some other function set, including a set of stepping (piecewise constant) functions, could be analogously used to define the BS for delayed dynamics. Then, all these functions could be used in the procedure of cross-validation for getting a more steady estimation of the BS.

II. BS IN ADDITIONAL BENCHMARK MODELS WITH TIME DELAYS
We have demonstrated the generalized BS approach with the circumstance of multiple equilibriums in the main text. In fact, it could be applicable to systems with stable periodic orbits and even systems with both kinds of stable steady states, which is illustrated by the following examples.
A. Time-delayed Van der Pol-Duffing oscillator First, we investigate the time-delayed Van der Pol-Duffing oscillator, where a linear term of delayed position feedback is induced [3]. The equation of motion with x denoting the position of an oscillator is where α and γ are damping coefficients, β is the rigidity coefficient, and ω 0 is the inherent frequency of the system. In the feedback term, x τ = x(t − τ) with τ denotes the time delay and A is the gain coefficient of feedback. Here, we set α = 8 and α = 15 for the expansions of the initial values of x andẋ, respectively, and denote the BS simply by S n B .
Coexistence of multiple periodic orbits and equilibriums of system (S1) can be found when the parameters are properly selected. In particular, we perform our BS computing approach with changing time delay τ in the interval [4,8]. As shown numerically in Fig  B for the periodic orbit of a larger amplitude for the three bases, which shows a high consensus.

B. Coupled Stuart-Landau oscillators with time delays
The Stuart-Landau oscillator also exhibits co-existence of multiple stable periodic orbits when delayed state feedback terms are induced, which is completely different from the system without delay [4]. We use the generalized BS method to investigate the variation of the basin stability with time delay for the co-existing stable periodic orbits. The Stuart-Landau system with delay feedback terms can be written as: where α, β, ω 0 are parameters and K is a gain coefficient. Here, we set α = 12 and denote by S n B = S (n,12)

III. STABILITY ANALYSIS FOR SELECTING THE RANGE OF COEFFICIENTS
In the design of the generalized BS for the delayed dynamics, the n-th order approximate BS, S (n,α) B , is clearly relevant to the parameter α. Different α may lead to different S (n,α) B . However, what we are mainly concerned about is the changing tendency S (n,α) B with respect to the essential parameters. Thus, it is adequate for us to choose sufficiently large α for calculating the BS. Moreover, for each coefficient a i in front of the basis functions, we still need to determine the range [γ 1 , γ 2 ] from which we can sample the values. Once all the sampled a i satisfy n i=1 a 2 i ≤ α 2 , the corresponding function is taken as an initial value for the delayed dynamics.
Consider a stable steady state, denoted by A. Then the corresponding sampling population becomes U = {Z 1 , Z 2 , · · · , Z N }, where Z i = 1 represents the i-th initial value in E (n,α) from which the trajectory of the delayed dynamics is attracted to A, and otherwise Z i = 0. We sample k functions from U and define an indicator function as follows: 1, when Z i is sampled, 0, when Z i is not sampled.
We intend to estimate the proportion of N 1 = N i=1 Z i in N, i.e., Z = N 1 N = p, which is equal to S (n,α) B . To this end, we need to use the following well-known Wald-Wolfowitz Theorem [5].
Wald-Wolfowitz Theorem. Suppose that {a N1 , · · · , a NN } and {x N1 , · · · , x NN } are two sequences of real numbers. Also suppose that these sequences satisfy for r = 3, 4 and larger N, where For each N, (X 1 , · · · , X N ) is supposed to be a random vector uniformly taking value from all the permutations of the sequence {x N1 , · · · , x NN } and denote by a Ni X i .
Then, we have Moreover, the probability for L N has the following convergence law: as N → ∞.
In order to use the theorem above for solving our problem, we choose It is easy to check that these settings satisfy the conditions required in the theorem above. In what follows, we compute E(L N ) and Var(L N ). On the one hand, On the other hand, where the last equality is due to as N → ∞. Here, N(0, 1) represents the standard normal distribution. Given a confidence level 1 − β, we have an interval estimation for p as follows: where u 1− β 2 is the 1 − β 2 -quantile of the distribution N(0, 1). Solving the inequality for p in (S3) yields: P(P L ≤ p ≤ P U ) ≈ 1 − β as N is sufficiently large, where The above argument manifests that for a sufficiently large N, the estimated interval [P L , P U ] for p = S (n,α)

B
with an increase of k, reduces to a fairly small interval in which p is located with a large probability.
As an illustrative example, we still use the time-delayed Van der Pol-Duffing oscillator (S1) with τ = 5. We adopt the Bernstein basis to estimate the interval of [P L , P U ] within a confidence level 1 − β =