Artur Rego-Costa

Florence Débarre

Luis-Miguel Chevin
email: luis-miguel.chevin@cefe.cnrs.fr

Chaos and the (un)predictability of evolution in a changing environment

Keywords: Eco-evolutionary dynamics, Adaptation to changing environments, Predictability, Repeatability, Chaotic dynamics

Among the factors that may reduce the predictability of evolution, chaos, characterized by a strong dependence on initial conditions, has received much less attention than randomness due to genetic drift or environmental stochasticity. It was recently shown that chaos in phenotypic evolution arises commonly under frequency-dependent selection caused by competitive interactions mediated by many traits. This result has been used to argue that chaos should often make evolutionary dynamics unpredictable. However, populations also evolve largely in response to external changing environments, and such environmental forcing is likely to influence the outcome of evolution in systems prone to chaos. We investigate how a changing environment causing oscillations of an optimal phenotype interacts with the internal dynamics of an eco-evolutionary system that would be chaotic in a constant environment. We show that strong environmental forcing can improve the predictability of evolution, by reducing the probability of chaos arising, and by dampening the magnitude of chaotic oscillations. In contrast, weak forcing can increase the probability of chaos, but it also causes evolutionary trajectories to track the environment more closely. Overall, our results indicate that, although chaos may occur in evolution, it does not necessarily undermine its predictability.

Introduction

The extent to which evolution is repeatable and predictable bears on the usefulness of evolutionary biology as a tool for a growing number of applied fields, including drug resistance management in pests and pathogens, or sustainable agriculture and harvesting under climate change. So far, the investigation of factors that may reduce the predictability of evolution has mostly focused on various sources of stochasticity (i.e. randomness), namely genetic drift, the contingency of mutations, and randomly fluctuating environments [START_REF] Crow | An Introduction to Population Genetics Theory[END_REF][START_REF] Lenormand | Stochasticity in evolution[END_REF][START_REF] Saether | The concept of fitness in fluctuating environments[END_REF]. A much less explored source of unpredictability in evolution is deterministic chaos (but see [START_REF] Hamilton | Sex versus Non-Sex versus Parasite[END_REF][START_REF] Altenberg | Chaos from Linear Frequency-Dependent Selection[END_REF][START_REF] Gavrilets | Intermittency and Transient Chaos from Simple Frequency-Dependent Selection[END_REF][START_REF] Abrams | Prey Adaptation as a Cause of Predator-Prey Cycles[END_REF]Dercole & Rinaldi, 2010;Dercole et al., 2010;[START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF], which occurs when the dynamics of a system, despite being completely predictable from the initial conditions, are strongly sensitive to them. Under chaotic dynamics, any measurement error, regardless how small, will be amplified over time, to the point that it becomes impossible to make accurate predictions beyond a certain timescale (Ott, 2002;[START_REF] Petchey | The ecological forecast horizon, and examples of its uses and determinants[END_REF].

It was recently demonstrated theoretically that evolutionary dynamics at the phenotypic level can become chaotic when natural selection results from between-individual interactions within a species [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF], i.e., even in the absence of any interspecific eco-evolutionary feedbacks, which are known to enhance ecological chaos [START_REF] Abrams | Prey Adaptation as a Cause of Predator-Prey Cycles[END_REF]Dercole & Rinaldi, 2010;Dercole et al., 2010). Specifically, evolutionary chaos arises in single-species models when (i) selection is frequency-dependent, such that the fitness of an individual depends on trait-mediated interactions with conspecifics; (ii) the fitness effects of these interactions are not simply a function of the phenotypic difference between the interactors (unlike typical models of within-and between-species interactions, e.g. [START_REF] Dieckmann | On the origin of species by sympatric speciation[END_REF][START_REF] Doebeli | Complexity and diversity[END_REF]; and (iii) the number of traits involved in these interactions (described as organismal complexity) is large [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF]. The authors concluded from this study that evolution is likely to be much less predictable than generally perceived. However, the theoretical demonstration that chaos is possible in a system is not sufficient to argue that this system is necessarily unpredictable, as we elaborate below.

First, the parameter values that lead to chaos may be rare in nature [START_REF] Hastings | Chaos in Ecology: Is Mother Nature a Strange Attractor?[END_REF][START_REF] Zimmer | Complex systems: Life After Chaos[END_REF]. For instance in ecology, chaos has long been known to be a possible outcome of even simple population dynamic models [START_REF] May | Simple mathematical models with very complicated dynamics[END_REF], but despite a few clear empirical demonstrations in the laboratory [START_REF] Benincà | Chaos in a long-term experiment with a plankton community[END_REF] and in the wild [START_REF] Benincà | Species fluctuations sustained by a cyclic succession at the edge of chaos[END_REF], most natural populations seem to have demographic parameters placing them below the "edge of chaos" [START_REF] Hastings | Chaos in Ecology: Is Mother Nature a Strange Attractor?[END_REF][START_REF] Ellner | Chaos in a Noisy World: New Methods and Evidence from Time-Series Analysis[END_REF][START_REF] Zimmer | Complex systems: Life After Chaos[END_REF][START_REF] Dercole | Coevolution of slow-fast populations: evolutionary sliding, evolutionary pseudo-equilibria and complex Red Queen dynamics[END_REF].

Second, and importantly with respect to evolution, a potentially chaotic system may still be predictable because it is subject to forcing by external factors with autonomous dynamics. In eco-evolutionary processes, such external forcing generally results from a changing environment that affects fitness components and their dependence on phenotypes, causing variation in natural selection. In fact, environmental variability affecting population growth [START_REF] Lande | Stochastic Population Dynamics in Ecology and Conservation[END_REF][START_REF] Ellner | Does rapid evolution matter? Measuring the rate of contemporary evolution and its impacts on ecological dynamics[END_REF][START_REF] Pelletier | Decomposing variation in population growth into contributions from environment and phenotypes in an age-structured population[END_REF] and natural selection [START_REF] Maccoll | The ecological causes of evolution[END_REF][START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF] is probably ubiquitous in natural populations, as documented notably by numerous examples of ecological and evolutionary responses to climate change (reviewed by [START_REF] Davis | Evolutionary Responses to Changing Climate[END_REF]Parmesan, 2006;[START_REF] Hoffmann | Climate change and evolutionary adaptation[END_REF]. When such external forcing is imposed, the predictability of evolution is likely to change, because (i) forcing can alter the probability that the system is indeed chaotic, for instance by suppression of chaos through synchronization to the external forcing (e.g. [START_REF] Pikovsky | Synchronization: A Universal Concept in Nonlinear Sciences[END_REF]; and, (ii) even if the dynamics remain chaotic, they may still be affected by the forcing in ways that make them largely predictable.

We investigate how a changing environment affecting phenotypic selection, modeled as a moving optimal phenotype (a classic approach, reviewed by [START_REF] Kopp | Rapid evolution of quantitative traits: Theoretical perspectives[END_REF], influences the predictability of evolutionary dynamics in a context where chaos is expected to make evolution highly unpredictable in the absence of environmental forcing. Focusing on a periodic environment, we ask how the amplitude and period of cycles influence (i) the probability that evolutionary trajectories are chaotic, and (ii) the degree to which those trajectories that are indeed chaotic track the optimal phenotype set by the environment, making them partly predictable. We show that a changing environment can dramatically alter the probability of evolutionary chaos in either direction, but that evolutionary tracking of the environmental forcing generally contributes to making evolutionary trajectories much more predictable than anticipated from a theory that ignores any role of the external environment. This suggests that the predictability of evolution is partly determined by a balance between the strength of intraspecific interactions and responses to environmental change.

Methods

Model

We consider a set of d phenotypic traits that evolve under both frequency-dependent and frequency-independent selection. Frequency-independent selection is assumed to be caused by stabilizing selection towards an optimal trait value, at which carrying capacity K is maximized. For instance, selection on beak size/shape in a bird, or mouth shape in a fish, may have a local optimum set by the available type of resources [START_REF] Martin | Multiple Fitness Peaks on the Adaptive Landscape Drive Adaptive Radiation in the Wild[END_REF][START_REF] Grant | 40 Years of Evolution: Darwin's Finches on Daphne Major Island[END_REF]. The frequency-dependent component of selection, on the other hand, emerges from trait-mediated ecological and social interactions between individuals within the species (either cooperative or competitive). Selection on a bird's beak morphology, for instance, depends not only on the available types of resources, but also on competition with conspecifics for these resources (as in [START_REF] Grant | 40 Years of Evolution: Darwin's Finches on Daphne Major Island[END_REF]. The intensity of this competition may depend on the beak size of competing individuals, but also on other traits of these interactors, such as their aggressiveness, territoriality, or degree of choosiness in food preference. This frequency-dependent component of selection, when it involves many traits, can lead to complex dynamics such as chaos or internally driven cycles, as shown by [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF].

Other details of the models, laid out in the Supporting Information, follow [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF] for ease of comparison. Notably, we use the same adaptive dynamics assumptions, whereby evolution is slower than population dynamics and relies on rare new mutations [START_REF] Dieckmann | The dynamical theory of coevolution: a derivation from stochastic ecological processes[END_REF][START_REF] Geritz | Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF][START_REF] Dercole | Analysis of evolutionary processes: the adaptive dynamics approach and its applications[END_REF], although our results are likely to apply also in a quantitative genetic context where evolutionary dynamics are much faster (see Discussion). Under these assumptions, the evolutionary dynamics of each phenotypic trait x i is (see Supporting Information for more details)

dx i dt = ∑ j = 1 d b i j x j + ∑ j, k = 1 d a i jk x j x k -x i -θ i t 3 . (1)
The first two terms in equation ( 1) represent frequency-dependent selection caused by phenotype-dependent interactions between individuals. The coefficients b ij and a ijk determine, respectively, the strength of first-and second-order selective interactions between traits. The latter occurs when the fitness of an individual with a given phenotype x i at trait i depends on the product of x i with the phenotypes at two traits j and k (including j = i and/or k = i), at least one of which is from its interactor (e.g. beak size of focal individual interacting with beak size and agressiveness of interactor). Importantly, these frequencydependent components of directional selection will be null if the interaction between individuals depends only on their phenotypic difference, in which case frequency-dependent selection would not lead to chaotic dynamics (as explained in the Supporting Information). Put differently, this means that a necessary (but not sufficient) condition for evolutionary chaos in this model is that intraspecific interactions do not depend solely on the resemblance (or difference) between the trait values of interactors, but also on their actual phenotype, for instance when individual that are highly social, agressive, or large, interact more overall.

The last term in equation (1) models stabilizing selection that causes the phenotype to evolve towards the optimum θ i for each trait i. In the original model [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF], the optimal phenotype was assumed to be constant and equal to zero for all traits. Here in contrast, we modeled the forcing effect of a changing environment by letting the optimum θ(t) change with time. We focused on a cycling environment causing the optimum to oscillate sinusoidally, which may represent, depending on the organism, oscillations in biotic (predators, parasites) or abiotic (e.g. meteorological) conditions on seasonal, multi-annual (e.g., El Niño oscillation), or geological time scales. We assume for simplicity that the optimal values for all traits respond to the same underlying environmental variable, such that they oscillate with the same period and phase (i.e., they are synchronized). The vector of optimal phenotypes for all traits can then be written as

θ t = A • sin 2π T t , ( 2
)
where T is the period, and A= (A 1 , A 2 , …, A d ) is a vector of amplitudes of oscillation for each of the d traits. This defines a single sine wave of amplitude ‖A‖ (the norm of vector A), and direction given by the unitary vector A/‖A‖. In our simulations, we focused for simplicity on the case where A i is the same for all traits, such that the optimum changes along a diagonal of the phenotype space. Note also that equation ( 2) implies that the optimum fluctuates around the phenotype for which the strength of selective interactions vanishes (set to the origin by definition, without loss of generality). Allowing for fluctuations to be centered on a different phenotype -or equivalently, including a 0 th order term in the interaction component of selection in equation ( 1) -would select for increasing interactions in all environments, thus artificially increasing the probability of chaos relative to [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF], where the optimum was set constant at the origin.

Simulations

We studied the evolution of the phenotypic vector x(t) containing the mean phenotype for each trait by numerically solving the dynamics of equation ( 1) given the initial phenotype x(0), period T and vector of amplitudes A (both characterizing the optimum θ), and set of interaction coefficients a ijk and b ij . In each simulation, the interaction coefficients were independently drawn from a normal distribution with mean 0 and standard deviation 1, and rescaled as

a i jk d and b i j d .
This rescaling, proposed by [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF], ensures that dynamics under different dimensionalities under a constant environment explore similar ranges of phenotypic values, between -1 and 1. It also prevents the unrealistically strong selection produced at high dimensionalities in the model with unscaled interaction coefficients: in effect, we keep the expected overall strength of selection constant but spread it across the d traits. Trajectories were run up to t = 1200 using the LSODA method, as implemented in the package deSolve in R [START_REF] Soetaert | Solving Differential Equations in R : Package deSolve[END_REF][START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF], with integration step dt = 0.1. Initial phenotypes were drawn from a multivariate normal distribution centered at zero such that, on average, the carrying capacity K(x(0), θ(0)) = 0.5 (Eq. ( S2)). This choice was made to keep the initial state of the system under biologically relevant degrees of adaptation.

Evolutionary dynamics were categorized based on their largest Lyapunov exponent λ, which measures the rate of exponential increase in the distance between two trajectories that start from very close initial conditions [START_REF] Sprott | Chaos and Time-Series Analysis[END_REF]. Dynamics that converge to an equilibrium phenotype have λ < 0, those that oscillate periodically (limit cycles) have λ = 0, and those that systematically diverge due to strong sensitivity to initial conditions (which defines chaos) have λ > 0. Here we used a local average Lyapunov exponent computed over a sliding window of 200 time units (see Supporting Information and Supplementary Fig. S1).

Proportion of transient chaos in a constant environment

In this model, trajectories that eventually reach fixed points or limit cycles may exhibit complex behaviors for long periods of time (Fig. 1A andB), during which they are indistinguishable from chaos. To understand the prevalence in the system of this so-called transient chaos [START_REF] Grebogi | Crises, sudden changes in chaotic attractors, and transient chaos[END_REF][START_REF] Gavrilets | Intermittency and Transient Chaos from Simple Frequency-Dependent Selection[END_REF][START_REF] Lai | Transient Chaos: Complex Dynamics on Finite-Time Scales[END_REF], we ran simulations with d ranging from 2 to 100 traits, under a constant optimum set at zero for all traits. Trajectories were identified as transiently chaotic if they switched from λ > 0 to either λ = 0 or λ < 0 (i.e. they went from chaos to a cycle or fixed equilibrium, respectively) by the end of the simulations (at time t = 1200; excluding the first 50 time steps of the simulations that correspond to a burn-in period).

Transience is expected to cause the proportion of chaos-like trajectories to decrease exponentially with time [START_REF] Yorke | Metastable chaos: The transition to sustained chaotic behavior in the Lorenz model[END_REF], even for high-dimensional systems [START_REF] Grebogi | Crises, sudden changes in chaotic attractors, and transient chaos[END_REF]. For each dimensionality, we thus estimated the proportion f(t) of trajectories behaving chaotically at each time step, and used this to estimate the asymptotic proportion of trajectories that remain chaotic over infinite time. This was done using the non-linear least-squares method (nls function in R's stats package [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] to fit a statistical model of the form f(t) = A exp(Bt) + C + ε, where A, B and C are the estimated variables (C being the asymptotic proportion of chaos), and ε is the residue (see Supplementary Fig. S2).

Predictability in a changing environment

To investigate the effect of a changing environment, we focused on a system of high dimensionality (d = 70), because this leads to a high probability of chaos in a constant environment (Supplementary Fig. S3 and [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF]. We used 100 sets of interaction coefficients b ij and a ijk and initial phenotypes x 0 . Each set of parameters was used for simulations in a constant environment and in changing environments (40 different combinations of amplitude ‖A‖ and period T of optimum oscillation, with ‖A‖ ∈ {1.30,2.33,3.73,5.04,5.99} and T ∈ {1.5,2,3,5,10,20,50,100}.) The amplitudes were chosen such that the smallest carrying capacity that a phenotype centered at the origin would experience (i.e. when the optimal phenotype was at the peak of its oscillation) was {0.99,0.9,0.5,0.1,0.01}.

For each regime of environmental change, we investigated the extent to which chaotic trajectories (i.e. trajectories with a final Lyapunov exponent λ > 0) track the moving optimum. For this, the time-series of phenotypic values were regressed on the oscillating optimal phenotype. We focused on the projection x̂ of the multivariate phenotypes along the direction of oscillations of the optimum in the phenotypic space, which we regressed on a similar projection θ ̂ for the optimum (which is simply the norm of θ). Additionally, because the evolving phenotype systematically lags behind the moving optimum in such a system [START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF], we maximized the R 2 of the regression by shifting forward the time-series of phenotype relative to that of the optimum (as shown in Supplementary Fig. S4). This maximum R 2 quantifies the proportion of variance in the evolutionary dynamics explained by (lagged) tracking of the optimum, so it is a measure of the predictability of evolution conditional on knowledge of the environment. This tracking component of evolution necessarily has the same period as the optimum (Supplementary Fig. S4), so approximating it as a sinusoidal function for simplicity, the regression slope of the phenotype x̂ on the (resynchronized) optimum θ ̂ is just the ratio of amplitudes of their cycles (as we confirmed in our simulations, see Supplementary Fig. S5).

We also looked for signatures of tracking of the changing environment in the evolutionary time-series x̂ using spectral analysis (fast Fourier transform method, as implemented in R's stats package, [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. This technique treats time-series as superpositions of sine and cosine waves of different frequencies [START_REF] Shumway | Time Series Analysis and Its Applications -With R Examples[END_REF], and allows estimation of the amplitudes associated to each frequency of oscillation (1/T). More precisely, the spectral density as computed by the spectrum() function in R is half the squared amplitude of waves of the corresponding frequency of oscillation.

Results

Transient evolutionary chaos is common

Chaos in phenotypic evolution was previously shown to be a common outcome in a constant environment, when interactions with conspecifics cause frequency-dependent selection as modeled in equation ( 1) and the interactions are mediated by many traits [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF][START_REF] Ispolatov | Chaos in high-dimensional dissipative dynamical systems[END_REF]. However, some of the apparently chaotic trajectories are actually transient (Fig. 1A andB). The fraction of trajectories that exhibit such transient chaos strongly depends on organismal complexity d (Fig. 1C), being highest for d between 40 and 70 (approximately), where the proportion of trajectories that are chaotic is intermediate (Fig. S2, and [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF][START_REF] Ispolatov | Chaos in high-dimensional dissipative dynamical systems[END_REF]. The statistical models fitted to the proportion of chaotic trajectories in time predict that virtually all transient trajectories should eventually transition to fixed points or periodic cycles by t = 3200 (Fig. 1C). However before this transition, the phenotype may undergo very complex, random-like dynamics, with no indication that they will eventually transition to a state that, once established, is quite predictable (see Supplementary Fig. S3).

A changing environment may either increase or decrease the probability of chaos

We next investigated the effect of a changing environment. We tracked the proportion of trajectories that were chaotic (up to time t = 1200) under sinusoidal cycles in the optimal phenotype, with varying periods and amplitudes. The proportion of chaotic trajectories was strongly influenced by the regime of optimum oscillation (Fig. 2A). Long periods of oscillation increased the probability of chaos relative to a constant environment, for the same organismal complexity (d = 70). For long enough periods, essentially all trajectories became chaotic at dimensionality d = 70, so the chaos-enhancing effect depended little on the amplitude of oscillations. However at smaller dimensionality (d = 40), the probability of chaos was maximized for a period that depended on the amplitude of oscillations, decreasing for very large periods when the amplitude was large (Supplementary Fig. S6).

In contrast to this chaos-enhancing effect, a combination of large amplitudes and short periods of optimum oscillation led to a sharp decrease in the proportion of chaotic trajectories (Fig. 2A). This decrease at time t = 1200 is not simply caused by an earlier transition out of chaos by transient trajectories: the observed proportion of chaotic trajectories can be substantially lower than that projected in infinite time for a constant optimum condition.

Part of the effect of the changing environment on the probability of chaos is explained by the absolute speed of the moving optimum. A simple metric of speed is the norm of the derivative of equation ( 2) relative to time, whose average over a period is equal to 4‖A‖/T. We use this average value as our measure of optimum speed. A slow optimum (speed below 5 for d = 70, Fig. 2B) increases the probability of chaos relative to a constant environment, while a fast optimum (speed above 5 for d = 70, Fig. 2B) decreases the probability of chaos.

Furthermore, conditions that lead to a reduced probability of chaos relative to a constant environment (frequent oscillations with broad amplitude) also correspond to strong environmental forcing on the evolutionary dynamics, causing substantial maladaptation and directional selection even in a context without frequency dependence [START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF], as shown in Supplementary Fig. S7.

Strong environmental forcing renders chaotic trajectories more predictable

The predictability of evolution is not entirely captured by the probability of chaos, even in a deterministic system as modeled here, because chaotic evolutionary trajectories need not be entirely unpredictable. In fact these trajectories, despite looking erratic, are constrained to remain near the optimal phenotype imposed by stabilizing selection [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF]. When a changing external environment causes movements of the optimum, evolution necessarily tracks this moving optimum to some extent, even for chaotic trajectories, as illustrated in Fig. 3A andB. The predictability of evolution for those trajectories that are chaotic thus depends on how the internally driven dynamics due to intraspecific interactions interplay with forcing by the external environment.

To investigate how tracking of the environment affects the predictability of evolution, we regressed the time series of the evolving phenotype on that of the optimum, for trajectories that are chaotic in a changing environment. Such regressions are similar to phenotypeenvironment associations, as are commonly estimated empirically from time series of phenotypes and environments (review in e.g. [START_REF] Michel | Evolution of phenotype-environment associations by genetic responses to selection and phenotypic plasticity in a temporally autocorrelated environment[END_REF], and similarly across space (e.g. [START_REF] Phillimore | Differences in spawning date between populations of common frog reveal local adaptation[END_REF]. Note that here, the phenotype first needs to be resynchronized with the environment, to correct for the adaptational lag (Supplementary Fig. S4, [START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF]. The R 2 of the regression of phenotype on the environment is a measure of the predictability of evolutionary trajectories conditional on knowledge of the environment, since it measures the proportion of the total temporal variance in phenotype that is explainable by tracking of the optimum. The predictability of chaotic evolutionary dynamics increases with larger amplitudes and longer periods of optimum oscillation (Fig. 4A; and see Supplementary Fig. S8 for regressions on optimum speed). Indeed, (i) evolutionary trajectories track long-period oscillating optima more closely than they do short-period ones [START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF], and (ii) if this optimum undergoes ample fluctuations, so will the phenotype, such that the tracking component of evolution will be substantial, as illustrated in Fig. 3B.

The exact same pattern is found for the repeatability of evolutionary trajectories among replicates, as measured by the proportion of the total variance in evolutionary trajectories attributable to temporal variation in the mean trajectory across replicates (Supplementary Fig. S9). The rationale for this measure is that the mean trajectory captures the tracking component of evolutionary change, while any additional variance around this trajectory is contributed by the internal dynamics caused by frequency-dependent selection. Repeatability thus has a very similar meaning to predictability, but unlike the latter it does not require information about the environment. Note also that repeatability is a measure of crosspredictability of evolution across replicates, and hence a metric of parallel phenotypic evolution, as commonly reported in the laboratory and in the field [START_REF] Lenormand | Parallel evolution: what does it (not) tell us and why is it (still) interesting?[END_REF].

The amplitude of the tracking component of phenotypic evolution (scaled to the amplitude of the optimum oscillation) is quantified by the regression slope of phenotype on optimum. This relative amplitude increases with the period of cycles in the optimum, tending towards 1 for long-period cycles of any amplitude (Fig. 4B). In contrast for shorter periods of optimum oscillation (left half of Fig. 4B), evolutionary tracking of the environment is more moderate, and less efficient for larger amplitudes.

The chaotic, non-predictable component of the evolutionary dynamics is captured by the residuals of the regression of the phenotype on the optimum. The variance of residuals was lower under combinations of large amplitudes and small periods of environmental change (Fig. 4C), corresponding to fast optima. Therefore, strong environmental forcing further increases the predictability of chaotic trajectories by dampening the extent to which they oscillate chaotically around the optimum. This effect can be seen in Fig. 3A, where the trajectory departs less from the optimum (due to chaos) than the one in Fig. 3B, with same amplitude but larger period.

Finally, we used spectral analysis, a general signal processing approach, to investigate whether environmental tracking leaves a signature in individual evolutionary time series.

Similarly to the repeatability analysis, this involves the investigation of the evolutionary dynamics without knowledge of the forcing environment. For essentially all regimes of environmental change that we investigated, the frequency of oscillation with highest spectral density in the time series of phenotypic evolution corresponded to that of the oscillating optimum, as illustrated for specific cases in Fig. 3C andD. The amplitude of the wave corresponding to this frequency of oscillation, as estimated from the spectral density, matched the amplitude estimated from the regression analysis (Fig. 4B) very well (Supplementary Fig. S5). This shows that even in a context where evolutionary dynamics can be largely chaotic, the strongest signal in the evolutionary time series is likely to be the predictable response to environmental change.

Discussion

When investigating factors that reduce predictability of evolution, evolutionary biologists have mostly focused on different sources of stochasticity [START_REF] Crow | An Introduction to Population Genetics Theory[END_REF][START_REF] Lenormand | Stochasticity in evolution[END_REF][START_REF] Saether | The concept of fitness in fluctuating environments[END_REF], while the possibility of chaos has received comparatively less attention. [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF] have recently shown that interactions mediated by many traits can produce chaos in phenotypic evolution. However, the implications for the predictability of evolution should not be overemphasized. First, the condition for chaos to occur is not simply frequency-dependent selection on many traits: it also requires that the strength of intraspecific interactions depends directly on the phenotypes of interactors, rather than only on their phenotypic match or distance [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF], while interaction of the latter type are more commonly used in models of evolutionary ecology. Second, we reveal that chaos can be transient in this model, such that initially chaotic dynamics might not be observed after some time. And third, forcing by a changing environment, a ubiquitous driver of eco-evolutionary dynamics (as exemplified by responses to climate change, [START_REF] Davis | Evolutionary Responses to Changing Climate[END_REF]Parmesan, 2006;[START_REF] Hoffmann | Climate change and evolutionary adaptation[END_REF], can modify the predictability of evolution in a system that may be chaotic in a constant environment. Below we discuss this latter point further.

Chaos vs. forcing

Strong environmental forcing (represented by a fast moving optimum) caused suppression of chaos relative to a constant environment in our model. In the general literature on chaos (outside of evolutionary biology), it has previously been demonstrated, both theoretically and experimentally, that chaotic dynamics can be suppressed with even slight forcing, but this usually occurs when the period of the forcing oscillation aligns closely to one of the specific resonance periods of the dynamics [START_REF] Lima | Suppression of chaos by resonant parametric perturbations[END_REF][START_REF] Fronzoni | Experimental evidence of suppression of chaos by resonant parametric perturbations[END_REF]. In contrast in our simulations, the proportion of chaotic trajectories was reduced over a broad range of short periods and large amplitudes of oscillation that jointly result in a fast moving optimum. In fact, most of the chaos suppression occurred for periods much shorter than the typical internal chaotic dynamics in a constant environment (Supplementary Fig. S10). One of the reasons for this discrepancy may be that we studied the effect of the moving optimum on a total of 250 different sets of interaction coefficients, instead of focusing on single defined chaotic system, as usually done in the physics literature. Another reason is that in our evolutionary model, the dynamical system is different (and in general more complex) than in models studied in physics, for instance as in [START_REF] Lima | Suppression of chaos by resonant parametric perturbations[END_REF] and [START_REF] Fronzoni | Experimental evidence of suppression of chaos by resonant parametric perturbations[END_REF].

In the opposite regime of slow optimum, the proportion of chaotic dynamics increased. Such an outcome is not uncommon in the general literature on chaos, as exemplified by the chaotic dynamics induced on a pendulum with friction by externally imposed sinusoidal impulses (forced damped pendulum, Ott, 2002). Such a phenomenon has also been shown in models of population ecology. [START_REF] Rinaldi | Multiple attractors, catastrophes and chaos in seasonally perturbed predator-prey communities[END_REF] explored seasonal oscillation of demographic parameters in a predator-prey model, and found that such forcing easily led to quasiperiodic and chaotic dynamics that did not occur in a constant environment. Similarly, [START_REF] Benincà | Species fluctuations sustained by a cyclic succession at the edge of chaos[END_REF] designed a model based on an empirical rocky shore ecosystem with non-transitive interactions of the rock-paper-scissor form, and showed that the model only led to the chaotic oscillations observed in nature when including the forcing effect of yearly seasonal temperature cycles on death rates. In our model, optima that move slowly cause repeated long excursions away from trait values where interactions are minimal (x = 0 in the model). This causes stronger interaction terms in equation (1) (and hence more chaos) relative to a constant optimum at x = 0, because chaos-enhancing interactions are stronger with larger trait values in this model [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF].

Generality and relevance to natural populations

Our results are difficult to compare quantitatively to empirical measurements, because we relied on the adaptive dynamics approach (to facilitate comparison with the original model in a constant environment by [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF], where time is not measured in units of time or in generations, but in terms of an implicit number of mutations fixed. In principle, this would suggest that the timescales of evolutionary and environmental change in our model are restricted to be slow, since they are limited by the time between independent fixation events. However, several authors [START_REF] Abrams | Evolutionarily unstable fitness maxima and stable fitness minima of continuous traits[END_REF][START_REF] Waxman | 20 Questions on Adaptive Dynamics[END_REF][START_REF] Débarre | Multidimensional (co)evolutionary stability[END_REF] have highlighted that the canonical equation of adaptive dynamics (on which eq. ( 1) is based) is almost identical to the equation for the response to selection in quantitative genetics [START_REF] Lande | Natural Selection and Random Genetic Drift in Phenotypic Evolution[END_REF], but with the latter operating on much shorter time scale. Furthermore, [START_REF] Ispolatov | Individual-based models for adaptive diversification in highdimensional phenotype spaces[END_REF] demonstrated that the chaotic attractors of the original [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF] paper can be closely matched by individual-based simulations that allow for substantial polymorphism, and thus for much faster evolutionary dynamics than under the classic adaptive dynamics assumptions (i.e., rare substitution events in otherwise monomorphic populations). This suggests that the relevant timescale of evolutionary and ecological change in our model can be much faster than assumed by adaptive dynamics, e.g. that of microbial evolutionary experiments or long-term field surveys, for which it is well known that evolutionary change can happen quickly (e.g. [START_REF] Hendry | The Pace of Modern Life: Measuring Rates of Contemporary Microevolution[END_REF][START_REF] Campbell-Staton | Winter storms drive rapid phenotypic, regulatory, and genomic shifts in the green anole lizard[END_REF]. Finally, understanding of the timescale of these chaotic oscillations will help elucidate the importance of transience, which should depend on the probability that an evolving transient population is observed before it has time to transition.

A crucial parameter of the model that is difficult to measure empirically is organismal complexity. The number of traits that can be measured is virtually infinite, but some may be highly correlated, or have negligible effects on fitness. From an evolutionary perspective, phenotypic complexity thus has to be defined with respect to its effects on fitness and selection. Previous theory has shown that complexity defined in a similar way as here (number of traits under stabilizing selection) has important impacts on the rate of adaptation [START_REF] Fisher | The Genetical Theory Of Natural Selection[END_REF][START_REF] Orr | Adaptation and the Cost of Complexity[END_REF], speciation and diversification [START_REF] Doebeli | Complexity and diversity[END_REF][START_REF] Chevin | Niche Dimensionality and the Genetics of Ecological Speciation[END_REF][START_REF] Débarre | Multidimensional (co)evolutionary stability[END_REF][START_REF] Svardal | Organismal complexity and the potential for evolutionary diversification[END_REF], or the drift load in a finite population [START_REF] Poon | Compensating for our Load of Mutations: Freezing the Meltdown of Small Populations[END_REF][START_REF] Tenaillon | Quantifying Organismal Complexity using a Population Genetic Approach[END_REF]. Some of these predictions have been used to attempt to estimate organismal complexity indirectly, through its emerging effects on fitness effects of mutations. Very different results have been obtained, with complexity ranging from very low (on the order of 1) to several orders of magnitude for the same organism, depending of the underlying model used [START_REF] Martin | A general multivariate extension of Fisher's geometrical model and the distribution of mutation fitness effects across species[END_REF][START_REF] Tenaillon | Quantifying Organismal Complexity using a Population Genetic Approach[END_REF]. In any case, our results under a changing environment should apply whenever the complexity of frequency-dependent selection caused by intraspecific interactions is high enough to generate chaotic dynamics in a constant environment [START_REF] Doebeli | Chaos and unpredictability in evolution[END_REF].

A possibility not explored in our study is that a changing environment alters the interaction strength between individuals with given phenotypes, and thus the frequency-dependent component of selection in equation (1). It is not entirely clear how to best model such an effect of the environment on interaction coefficients -while a moving optimum for the frequency-independent component of selection is well-established in models (reviewed by [START_REF] Kopp | Rapid evolution of quantitative traits: Theoretical perspectives[END_REF], and has some empirical support [START_REF] Chevin | Estimating the variation, autocorrelation, and environmental sensitivity of phenotypic selection[END_REF]. We can however anticipate that environments leading to larger absolute values of the interaction coefficients should increase the probability of chaos.

Interacting sources of unpredictability in evolution

We have focused for simplicity on effectively infinite populations in a fully deterministic environment, such that any unpredictability in evolution has to come from sensitivity to initial conditions characteristic of chaotic dynamics. This is perhaps a good approximation for some experimental studies with microbes in controlled environments, but more generally populations in the wild should also be exposed to another source of unpredictability: evolutionary stochasticity caused by genetic drift, the contingency of mutations, or a randomly changing environment [START_REF] Crow | An Introduction to Population Genetics Theory[END_REF][START_REF] Lenormand | Stochasticity in evolution[END_REF][START_REF] Saether | The concept of fitness in fluctuating environments[END_REF]. A more complete understanding of the predictability of evolution would therefore require combining stochasticity and chaos to investigate their possible interactions, as advocated previously for population dynamics [START_REF] Ellner | Chaos in a Noisy World: New Methods and Evidence from Time-Series Analysis[END_REF]. For instance, stochastic factors have been shown to either increase or decrease (depending on the chaotic system) the time that transients take to converge to their equilibrium states [START_REF] Lai | Transient Chaos: Complex Dynamics on Finite-Time Scales[END_REF].

While such an analysis is beyond our scope here, some preliminary statements can be made based on our results and those from the literature. Beyond just adding random variation among replicates (and thus directly reducing evolutionary predictability), stochasticity may interact with chaos in evolutionary dynamics, amplifying or reducing its importance. A stationary stochastic environment is a type of forcing that shares some similarities with deterministic cycles. Indeed, quantitative genetic models [START_REF] Lynch | Evolution and extinction in response to environmental change[END_REF][START_REF] Lande | The role of genetic variation in adaptation and population persistence in a changing environment[END_REF] have shown that increasing the stationary variance (respectively autocorrelation) of a stochastic environment has similar effects on the lag load (caused by phenotypic mismatches with the optimum) as increasing the amplitude (respectively period) of cycling environment. The results we report here could thus be used to guide interpretation about the outcomes from future models investigating the combined effects of stochasticity and chaos on the predictability of evolution. An interesting challenge for empirical research would be to establish the time scales at which chaos versus stochasticity dominate as sources of unpredictability in evolution. The linear regression of phenotype on the moving optimum gives insights about the predictability of chaotic evolutionary dynamics. (A) The fraction of the total temporal variation in the evolving phenotype explained by movements of the optimum, as captured by the R 2 of the regression, is higher when the optimum oscillates with larger amplitudes and longer periods. (B) The ratio of amplitudes between the tracking component of oscillations in the evolving phenotype and oscillations in the optimum, as estimated by slope of the regression, is smaller under short-period optimum oscillations. (C) Temporal variation in evolutionary dynamics introduced by chaos, beyond the variation attributable to tracking of the moving optimum, is captured by the variance of residuals in the regression. Conditions that allow for close tracking of the optimum in (B) also lead to dampened chaotic oscillations. We show the average (lines and points) and standard error (shading) over 100 simulations for each combination of amplitude and period of optimum oscillation.

Fig. 1 .

 1 Fig. 1.Transients in a constant environment. Transient trajectories present chaos-like behavior for some time, before transitioning to either (A) fixed equilibrium phenotypes or (B) periodic cycles. Two representative trajectories for a single trait are shown, simulated as described in the Methods, with d = 45 and a constant optimum (orange line). (C) The expected proportion of evolutionary trajectories (among all those for a given dimensionality d) that will eventually transition to non-chaotic dynamics, but would still be categorized as chaotic at a given time, are shown for different times. These proportions were estimated from a

Fig. 2 .

 2 Fig. 2.Probability of chaos in a changing environment. The proportion of chaotic trajectories at time t = 1200 in an oscillating environment is shown against (A) the period and (B) the average speed of optimum oscillation, for different values of amplitude (colors). For each condition of optimum oscillation (period and amplitude), we report the estimate (dot, line) and standard error (shading) of the proportion of chaotic trajectories, out of 100 replicated simulations that used sets of parameters also used in the constant-environment simulations (with dimensionality d = 70). The average speed of optimum oscillation is calculated from equation (2) as 4‖A‖/T. The observed proportion of chaos at t = 1200 (horizontal solid line), and the predicted proportion after all chaotic transients have transitioned (dashed line; calculated as shown in Supplementary Fig.S2) are also represented for a constant environment with d = 70. While long periods cause the proportion of chaotic trajectories to increase relative to a constant environment, short periods of large amplitudes (fast oscillations) cause this proportion to decrease, even below what would be predicted through more rapid transition of transiently chaotic trajectories.

Fig. 3 .

 3 Fig. 3.Chaotic evolutionary dynamics in changing environments. Chaotic evolutionary trajectories (black in A and B) combine internally driven chaos with external environmental forcing through tracking of the optimum. Chaotic oscillations can be either (A) longer or (B) shorter than the oscillations of the optimum (shown in orange), depending on the periods of the latter (3 and 50, respectively). (C, D) Spectral analysis of these same time-series of evolving phenotype (using all traits as described in the Methods) show a peak of spectral density (blue star) at the frequency of oscillation that corresponds to the oscillating optimum.

Fig

  Fig. 4.

Evolution. Author manuscript; available in PMC 2018 May 18.

Acknowledgements

We thank R. Gomulkiewicz, S. de Monte, and J. Pantel for insightful discussions, as well as C. Rueffler and four anonymous reviewers for useful feedback on the manuscript.

Funding

ARC is funded by an Erasmus Mundus Joint Master Degree scholarship from the Education, Audiovisual and Culture Executive Agency. FD is funded by a grant from Agence Nationale de la Recherche (ANR-14-ACHN-0003-01). LMC is funded by a grant from the European Research Council (ERC-2015-STG-678140-FluctEvol).

Conflict of interest:

The authors declare no conflict of interest.

Supplementary Material

Refer to Web version on PubMed Central for supplementary material.


Jared Miller
email: jarmiller@control.ee.ethz.ch

Matteo Tacchi
email: matteo.tacchi@gipsa-lab.fr

Didier Henrion
email: henrion@laas.fr

Mario Sznaier
email: msznaier@coe.neu.edu

Unsafe Probabilities and Risk Contours for Stochastic Processes using Convex Optimization

This paper proposes an algorithm to calculate the maximal probability of unsafety with respect to trajectories of a stochastic process and a hazard set. The unsafe probability estimation problem is cast as a primal-dual pair of infinite-dimensional linear programs in occupation measures and continuous functions. This convex relaxation is nonconservative (to the true probability of unsafety) under compactness and regularity conditions in dynamics. The continuous-function linear program is linked to existing probability-certifying barrier certificates of safety. Risk contours for initial conditions of the stochastic process may be generated by suitably modifying the objective of the continuous-function program, forming an interpretable and visual representation of stochastic safety for test initial conditions. All infinite-dimensional linear programs are truncated to finite dimension by the Moment-Sum-of-Squares hierarchy of semidefinite programs. Unsafe-probability estimation and risk contours are generated for example stochastic processes.

Introduction

This work performs risk analysis of stochastic processes by computing the probability that trajectories will reach an unsafe set X u . Trajectories evolve over a maximal time horizon of T in a state space X ⊆ R n , beginning in an initial set X 0 ⊆ X and possibly entering the unsafe set X u ⊆ X. We consider stochastic processes described by a generator L (e.g., a Stochastic Differential Equation (SDE)), for which X t is the time-dependent state probability distribution of the process starting from X 0 and terminating upon exit from X (with exit time distribution τ X = inf{t : X t ∈ ∂X}). The worst-case probability of unsafety for L-trajectories starting at a point x 0 ∈ X within a time horizon of [t 0 , T ] is:

P * (t 0 , x 0 ) = sup t * ∈[t0,T ] Prob X t * [x ∈ X u ] (1a)
x(t) follows L ∀t ∈ [t 0 , min(t * , τ X )] (1b)

x(0) = x 0 . (1c)
The stopping time min(t * , τ X ) in (1b) also ensures that trajectories x(t) remain in X for all relevant t. The worst-case probability of unsafety for L-trajectories over an initial set X 0 ⊆ X is:

P * (t 0 , X 0 ) = sup x0∈X0 P * (t 0 , x 0 ) (2a) = sup t * ∈[0,T ], x0∈X0 Prob X t * [x ∈ X u ] (2b)
x(t) follows L ∀t ∈ [t 0 , min(t * , τ X )] (2c)

x(0) = x 0 . (2d)
2 Preliminaries

Notation

The subset of natural numbers between a and b is a..b ⊂ N. The minimum of two quantities will be denoted as a ∧ b = min(a, b). The set of polynomials with real-valued coefficients in an indeterminate x ∈ R n is R[x]. Every polynomial c ∈ R[x] may be uniquely described by a finite sum over multi-indices α ∈ N n by c(x) = α c α x α . The degree of a polynomial deg c is equal to the largest exponent-sum n i=1 α i such that c α ̸ = 0. The set of polynomials in x with degree at most 2k is R[x] ≤2k .

The set of continuous functions over a set S is C(S). Its subcone of nonnegative continuous functions is C + (S). Given a product set S × H, the set of functions that are once-continuously differentiable in the first variable s and are twice-continuously differentiable in the second variable h is C 1,2 (S × H). The set of signed Borel measures supported in S is M(S), and its subset of nonnegative Borel measures supported in S is M + (S). The support of a measure µ ∈ M + (R n ) is the smallest closed set of points s such that every open neighborhood N ϵ (s) has positive measure µ(N ϵ (s)) > 0. The sets C + (S) and M + (S) are in topological duality when S is compact, and they admit a duality pairing ⟨•, •⟩ by Lebesgue integration: ∀f ∈ C + (S), µ ∈ M + (S) : ⟨f, µ⟩ = S f (s)dµ(s). This duality pairing will be extended into a bilinear pairing between C(S) and M + (S). The µ-measure of a set A ⊆ S may also be expressed as the pairing of µ with the indicator function

I A : µ(A) = ⟨I A , µ⟩ = inf{⟨w, µ⟩ | w ∈ C(S), w ≥ I A }.
The mass of µ is µ(S) = ⟨1, µ⟩, and µ is a probability distribution if this mass is 1. A vital example of a probability distribution is the Dirac delta δ s ′ supported only at s ′ , obeying the point-evaluation pairing rule ∀f ∈ C(S) : ⟨f,

δ s ′ ⟩ = f (s ′ ). The Lebesgue (volume) distribution over a space S is λ S with ∀f ∈ C(S) : ⟨f, λ S ⟩ = S f (s)ds. The product measure between µ ∈ M + (S 1 ) and ν ∈ M + (S 2 ) is the unique measure µ ⊗ ν satisfying ∀A 1 × A 2 ⊆ S 1 × S 2 : µ ⊗ ν(A 1 × A 2 ) = µ(A 1 )ν(A 2 ).

Stochastic Processes

A stochastic process is a time-indexed set of random variables {X t } that are related together through system dynamics [START_REF] Stroock | Multidimensional Diffusion Processes[END_REF] (pushforward of an initial distribution along flow maps). Properties of the stochastic process will be analyzed in terms of test functions in an appropriate set C (such as C = C([0, T ]×X). The expectation of a test function v(s, x) at the time t according to the distribution

X t is E[v(t, x) | X t ].
Letting ∆t > 0 be a time step, the generator L of a stochastic process satisfies (for all appropriate test functions v ∈ domL)

L ∆t v(t, x) = lim ∆t ′ →∆t E[v(t + ∆t ′ , x) | X t+∆t ′ ] -v(t, x) ∆t ′ . (3)
We will express the domain of L as C = domL, such that C is a subset of the preimage of continuous functions under L. The generator for a discrete-time stochastic process is L ∆t with ∆t > 0, which is defined w.r.t. the test function class C = C([0, T ] × X). For a discrete-time law following x t+∆t = f (t, x t , λ t ) in which the time-varying parameter λ t ∈ Λ has a probability distribution of ξ(λ t ), the associated generator satisfies

L ∆t v(t, x) = (∆t -1 ) Λ v (t + ∆t, f (t, x, λ)) dξ(λ) -v(t, x). (4)
The generator for a continuous-time (Feller) stochastic process is L ∆t=0 for the class C = C 1,2 ([0, T ]×X) [START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF].

The dynamical behavior of an Itô SDE has a description

dx = f (t, x)dt + g(t, x)dW, (5)
in which W is the Wiener process. The generator of the SDE in ( 5) is

L 0 v(t, x) =∂ t v(t, x) + f (t, x) • ∇ x v(t, x) + g(t, x) T ∇ 2 xx v(t, x)g(t, x)/2. (6)
A sequence of random variables [START_REF] Stroock | Multidimensional Diffusion Processes[END_REF] (in the sense of generalized conditional expectation). The generator L solves a martingale problem for all possible time steps s > 0 and test functions v ∈ C [START_REF] Williams | Probability with Martingales[END_REF]. The martingale problem in discrete-time (with s ∈ N) is

{Y t } is a martingale if E[Y t+δt | {Y t ′ } t ′ ≤t ] = Y t
0 = E[v(t + s, x) | X t+s ] -E[v(t, x) | X t ] (7)
-

t+s s ′ =t E[L ∆t v(t, x) | X s ′ ], and in continuous time (with s ∈ R ≥0 ) is 0 = E[v(t + s, x) | X t+s ] -E[v(t, x) | X t ] (8)
-

t+s s ′ =t E[L 0 v(t, x) | X s ′ ].
For the rest of this paper, we will refer to the generator L interchangeably with its stochastic process.

The martingale problem can be described using the theory of occupation measures. Given an initial condition (point) x 0 ∈ X 0 ⊆ X and a time t ∈ [0, T ], we refer to x(t | x 0 ) as the trajectory of {X t } with generator L starting at the initial distribution δ x=x0 . The occupation measure µ ∈ M + ([0, T ] × X) and terminal measure µ τ ∈ M + ([0, T ] × X) of the stochastic process {X t } from a starting time t 0 up to a stopping time t * satisfies ∀A ⊆ [t 0 , T ], B ⊆ X:

µ(A × B) = X0 t * t=t0 I A×B (t, x(t | x 0 )) dt dµ 0 (x 0 ) (9a) µ τ (A × B) = X0 I A×B (t * , x(t * | x 0 )) dµ 0 (x 0 ). ( 9b
)
The measures in (9) may also be defined with respect to a probability distribution of stopping times 2 obeys the martingale relations [START_REF] Abate | Probabilistic reachability and safety for controlled discrete time stochastic hybrid systems[END_REF] or (8) of the stochastic process with generator

t * ∈ M + ([0, T ]). A tuple of measures (µ 0 , µ, µ τ ) ∈ M + (X 0 ) × M + ([0, T ] × X)
L ∀v ∈ C : ⟨v, µ τ ⟩ = ⟨v(t 0 , •), µ 0 ⟩ + ⟨Lv, µ⟩. (10)
Relation ( 10) is called Liouville equation for Ordinary Differential Equations (ODEs) and is referred to as Dynkin's equation in the context of SDEs. We will also equivalently express [START_REF] Prajna | Stochastic Safety Verification Using Barrier Certificates[END_REF] in a shorthand form by

µ τ = δ t0 ⊗ µ 0 + L † µ. ( 11
)

Sum-of-Squares Methods

Let c ∈ R[x] ≤2k be a polynomial. The polynomial c is nonnegative if ∀x ∈ R n : c(x) ≥ 0. Verification of polynomial nonnegativity is generically an NP-hard task. One tractable method to construct nonnegative polynomials is through SOS representations. The polynomial c is SOS if there exists N polynomials q j ∈ R[x] such that c(x) = N j=1 q 2 j (x). The cone of SOS polynomials is Σ[x], and its subset of SOS polynomials with degree ≤ 2k is Σ[x] ≤2k . The set of SOS polynomials equals the cone of nonnegative polynomials only when n = 1, 2k = 2, or (n = 2, 2k = 4) [START_REF] Hilbert | Über die Darstellung definiter Formen als Summe von Formenquadraten[END_REF]. SOS polynomials become vanishingly small in the set of nonnegative polynomials as n and d rise [START_REF] Blekherman | There are Significantly More Nonnegative Polynomials than Sums of Squares[END_REF].

SOS representations up to fixed degree can be conducted by formulating SDPs. To each polynomial c ∈ R[x] ≤2k , there exists a (typically nonunique) associated polynomial vector [START_REF] Choi | Sums of squares of real polynomials[END_REF].

v k (x) ∈ (R[x] ≤d ) s and a symmetric Gram matrix Q ∈ S s such that c(x) = v k (x) T Qv k (x). One such choice of vectors v k (x) is the set of monomials in x between degrees 0..k, in which case s = n+k k . A polynomial c is SOS if its Gram matrix Q is also Positive Semidefinite (PSD)
A Basic Semialgebraic (BSA) set K is a set described by a finite number of bounded-degree polynomial inequality (N g ) and equality constraints (N h ):

K = {x ∈ R n | g i (x) ≥ 0, h j (x) = 0, ∀i ∈ 1..N g , j ∈ 1..N h }. ( 12
)
The class of Weighted Sum of Squares (WSOS) polynomials over K in [START_REF] Salamati | Data-driven safety verification of stochastic systems via barrier certificates[END_REF] is the set Σ[K] of polynomials c that admit a description of

c(x) =σ 0 (x) + i σ i (x)g i (x) + j ϕ j (x)h j (x) (13a) σ 0 , σ i ∈ Σ[x], ϕ j ∈ R[x]. (13b)
The set K is compact if there exists a finite

R > 0 such that {x | ∥x∥ 2 2 ≤ R} ⊇ K. The set is additionally Archimedean if R -∥x∥ 2 2 ∈ Σ[K].
Not every compact set is Archimedean [START_REF] Cimprič | Closures of quadratic modules[END_REF], but if an R verifying compactness of K is known, then the constraint R -∥x∥ 2 2 ≥ 0 can be added to the description of K to render K Archimedian. Every compact polytope and ellipsoid is also Archimedean [START_REF] Lasserre | Moments, Positive Polynomials And Their Applications[END_REF].

The Putinar Positivestellensatz states that every positive polynomial over K is also a member of Σ[K] when K is Archimedean [START_REF] Putinar | Positive Polynomials on Compact Semi-algebraic Sets[END_REF]. However, the polynomial degrees of the multipliers (σ, ϕ) needed to certify nonnegativity of such a positive c could be exponential in n and k [START_REF] Nie | On the complexity of Putinar's Positivstellensatz[END_REF].

Unsafe Linear Programs

This section will present LPs to analyze safety of a stochastic process with generator L.

Assumptions

We will posit the following assumptions throughout this paper: A1 The state-sets X 0 , X u , X are all compact.

A2 The time interval [t 0 , T ] is compact.

A3 Trajectories stop upon the first exit from X (τ X ∧ T ).

A4 The test function set

C = dom(L) satisfies C ⊆ C([t 0 , T ] × X) with 1 ∈ C and L1 = 0.
A5 The set C separates points and is multiplicatively closed.

A6 There exists a countable set {v

k } ∈ C such that ∀v ∈ C : (v, Lv) is contained in the bounded pointwise closure of the linear span of {(v k , Lv k )}.
Assumptions A4-A6 are the same as in Condition 1 of [START_REF] Cho | Linear Programming Formulation for Optimal Stopping Problems[END_REF]. We note that generators for Markov, SDE, and Lévy processes evolving in a compact set satisfy A1-A6.

Measure Program

We will pose a convex but infinite-dimensional LP in order to upper-bound program [START_REF] Esfahani | The stochastic reach-avoid problem and set characterization for diffusions[END_REF]. Such an LP will be posed in terms of a relaxed occupation measure (µ 0 , µ τ , µ).

Theorem 3.1. The following LP will upper-bound (2) with p * (t 0 , X 0 ) ≥ P * (t 0 , X 0 ) for each X 0 ⊆ X p * (t 0 , X 0 ) = sup ⟨I Xu , µ τ ⟩ (14a)

µ τ = δ t0 ⊗ µ 0 + L † µ (14b) ⟨1, µ 0 ⟩ = 1 (14c) µ 0 ∈ M + (X 0 ) (14d) µ, µ τ ∈ M + ([t 0 , T ] × X). (14e)
Proof. Let x 0 ∈ X 0 be an initial condition and t * ∈ [t 0 , T ] be a stopping time. A set of measures (µ 0 , µ, µ τ ) may be constructed to fulfill constraints (14b)-(14e). To begin, the initial measure µ 0 may be chosen as the Dirac delta δ x=x0 . The measure µ can be set to the occupation measure of the stochastic trajectory x(t | x 0 ) between the initial time of t 0 and the stopping time of t * ∧ τ X . The terminal measure µ τ (t, x) may then be picked as the the time-state distribution at time t * (in which trajectories are stopped according to t * ∧ τ X ).

Because there exists a mapping from every feasible point (t * , x 0 ) of (2) to a set of measures in (16b)-(16f), it holds that p * τ (t 0 , X 0 ) ≥ P * (t 0 , X 0 ).

The function I Xu (x) in the objective (14a) is a discontinuous function of x. We can reformulate ( 14) as an LP where all terms are continuous in (t, x) (under A1-A6) through the use of the following restriction lemma: Lemma 3.2 (Theorem 3.1 of [START_REF] Henrion | Approximate Volume and Integration for Basic Semialgebraic Sets[END_REF]). Let S 1 , S 2 be spaces with S 1 ⊆ S 2 , and let ξ ∈ M + (S 2 ) be a measure. The restriction of ξ to S 1 is uniquely given by the measure ν * solving the following optimization problem

Q * = sup ν,ν ⟨1, ν⟩ (15a)
ξ = ν + ν (15b) ν ∈ M + (S 1 ) (15c) ν ∈ M + (S 2 ). ( 15d
)
The objective Q * will satisfy Q * = ⟨1, ν * ⟩ = ⟨I S1 , ξ⟩. The measure ν * will be supported on closure(S 2 \ S 1 ).

We will define a peak measure µ p ∈ M + ([t 0 , T ]×X u ) and a complement measure µ c ∈ M + ([t 0 , T ]×X) in order to split the terminal measure µ τ as µ τ = µ p + µ c . This decomposition (through Lemma 3.2) will allow for the definition of an LP with continuous coefficients in (t, x). This continuity provides strong duality (under A1-A6) and convergence of the moment-SOS hierarchy to P * (t 0 , X 0 ) from (2) (under A1-A6 and assumptions on polynomial structure). The LP in which µ τ is decomposed as µ p + µ c is:

Theorem 3.3.
The following LP has the same objective value p * (t 0 , X 0 ) as in (2):

p * (t 0 , X 0 ) = sup ⟨1, µ p ⟩ (16a) µ p + µ c = δ t0 ⊗ µ 0 + L † µ (16b) ⟨1, µ 0 ⟩ = 1 (16c) µ 0 ∈ M + (X 0 ) (16d) µ, µ c ∈ M + ([t 0 , T ] × X) (16e)
µ p ∈ M + ([t 0 , T ] × X u ). (16f)
Proof. 1) Let (µ 0 , µ p , µ c , µ) be a feasible point to (16b)-(16f). Then (µ 0 , µ p + µ c , µ) is a feasible point to (14b)-(14e).

2) Let (µ 0 , µ τ , µ) be a feasible point to the constraints (14b)-(14e). Unique measures µ p , µ c may be defined according to the restriction Lemma 3.2 with S 1 = [t 0 , T ] × X u and S 2 = [t 0 , T ] × X by solving sup ⟨1, µ τ ⟩ (17a)

µ τ = µ p + µ c (17b)
µ c ∈ M + ([t 0 , T ] × X) (17c)
µ p ∈ M + ([t 0 , T ] × X u ), (17d)
with ⟨I Xu , µ τ ⟩ = ⟨1, µ p ⟩.

There exists a mapping relating each feasible point of ( 14) to ( 16) and vice versa while keeping the same objective, which proves that the two objectives are equal. Theorem 3.4. Under assumptions A1-A6, the objectives of ( 16) and (1b) are equal (p * (t 0 , X 0 ) = P * (t 0 , X 0 )).

Proof. By Theorem 3.1 of [START_REF] Cho | Linear Programming Formulation for Optimal Stopping Problems[END_REF] (and under Assumptions A1-A6), every tuple (µ 0 , µ, µ τ ) with µ 0 ∈ M + ([t 0 , T ]× X) (with ⟨1, µ 0 ⟩ = 1) and µ, µ τ ∈ M + ([t 0 , T ] × X) satisfying a Martingale relation [START_REF] Clark | Control barrier functions for stochastic systems[END_REF] is supported on the graph of a stochastic process (1b).

The tuple (µ 0 , µ, µ p + µ c ) from the Martingale relation ( 11) is supported on the graph of a stochastic process. The measures µ p and µ c can be separated according to Lemma 3.2. This concludes the no-relaxationgap proof.

Function Program

Theorem 3.5. An LP in terms of an auxiliary function v that forms a weak dual (p * (t 0 , X 0 ) ≥ d * (t 0 , X 0 )) to ( 16) is

d * (t 0 , X 0 ) = inf γ∈R γ (18a) γ ≥ v(t 0 , x) ∀x ∈ X 0 (18b) Lv(t, x) ≤ 0 ∀(t, x) ∈ [t 0 , T ] × X (18c) v(t, x) ≥ 0 ∀(t, x) ∈ [t 0 , T ] × X (18d) v(t, x) ≥ 1 ∀(t, x) ∈ [t 0 , T ] × X u (18e) v ∈ C. ( 18f
)
Strong duality (p * (t 0 , X 0 ) = d * (t 0 , X 0 )) will hold under assumptions A1-A4.

Proof. See Appendix A.

Remark 1. Constraints (18b)-(18f) have the same form as a stochastic (time-dependent) barrier function at a fixed probability level γ (from [START_REF] Prajna | A Framework for Worst-Case and Stochastic Safety Verification Using Barrier Certificates[END_REF]). Our work involves application of γ as an optimization variable, as well as a proof of no relaxation gap in Theorem 3.4 under A1-A6.

Risk Contour Linear Programs

The LPs in ( 16) and (18) determine the worst-case probability of unsafety among stochastic trajectories starting from an initial set X 0 . The optimal value of ( 16) and ( 18) is therefore equal to 1 when X 0 = X and X u is nonempty (because then ∅ ̸ = X u ⊂ X 0 and there exists trajectories starting in the unsafe set). This section will form a continuous-function LP similar to [START_REF] Ahbe | Region of attraction analysis of nonlinear stochastic systems using polynomial chaos expansion[END_REF], whose function values evaluated at a given initial point X 0 form a upper-bound on the probability of unsafety when starting at x 0 . These upper-bounds will be proven to converge in an L 1 sense over X to the true risk function P * . Such upper-bounds on unsafety allow for interpretation of risk when performing path planning tasks. Let µ 0 ∈ M + (X) be a fixed probability distribution (⟨1, µ 0 ⟩ = 1). The µ 0 -averaged probability of unsafety is

M * (t 0 , µ 0 ) = X P * (t 0 , x 0 ) dµ 0 (x 0 ). ( 19
)
One specific case we will focus on in the examples is where µ 0 is the the uniform distribution over X (µ 0 = λ X /vol(X)).

Measure Program

Corollary 1. An LP to upper-bound (2) with m * (t 0 , µ 0 ) ≥ M * (t 0 , µ 0 ) in general, and m

* (t 0 , µ 0 ) = M * (t 0 , µ 0 ) under A1-A6, is m * (t 0 , µ 0 ) = sup ⟨1, µ p ⟩ (20a) µ p + µ c = δ t0 ⊗ µ 0 + L † µ (20b) µ, µ c ∈ M + ([t 0 , T ] × X) (20c)
µ p ∈ M + ([t 0 , T ] × X u ). ( 20d
)
Proof. These optimality properties follow from arguments used in Theorem 3.3 (upper-bound by a measure construction starting from the distribution µ 0 ) and Theorem 3.4 (restriction measures). The measures (µ, µ p + µ c ) are supported on the graph of the stochastic process L starting from µ 0 for every feasible (µ, µ p , µ c ).

Function Program

For given v ∈ C and µ 0 ∈ M + (X), let v(t 0 , •) be shorthand notation for the function x → v(t 0 , x) such that ⟨v(t 0 , •), µ 0 ⟩ = ⟨v(t 0 , x), µ 0 (x)⟩. We will modify terms (18a)-(18b) in order to form an LP penalized over the whole space X with an objective integrating against the weighting distribution µ 0 :

J * (t 0 , µ 0 ) = inf⟨v(t 0 , •), µ 0 ⟩ (21a) Lv(t, x) ≤ 0 ∀(t, x) ∈ [t 0 , T ] × X (21b) v(t, x) ≥ 0 ∀(t, x) ∈ [t 0 , T ] × X (21c) v(t, x) ≥ 1 ∀(t, x) ∈ [t 0 , T ] × X u (21d) v ∈ C. (21e)
Lemma 4.1. The objective (21a) from ( 21) is upper-bounded by J * (t 0 , µ 0 ) ≤ 1.

Proof. The function v(t, x) = 1 is feasible for all constraints (21b)-(21e). Given that µ 0 is a probability measure, the objective in (21a) will satisfy J * ≤ ⟨1, µ 0 ⟩ = 1.

We now show that any solution of (21b)-(21e) returns an upper-bound of the probability of unsafety.

Theorem 4.2. Let v be feasible for (21b)-(21e) and let assumptions A1-A6 hold. Then

v(t 0 , x 0 ) ≥ P (t 0 , x 0 ) ( 22
)
for every initial condition x 0 ∈ X.

Proof. See Appendix B.

Remark 2. Level sets of any v(t, x) solving (21b)-(21e) can therefore be used to upper-bound the risk of stochastic execution for L when starting at any (t,

x 0 ) ∈ [t 0 , T ] × X.
The notation p(t 0 , x 0 ) will be subsequently used to denote the function obtained by solving ( 16) as x 0 → p(t 0 , X 0 = {x 0 }).

(t 0 , •), µ 0 ⟩ → J * (t 0 , µ 0 ) (under A1-A6). Then v k converges in L 1 (µ 0 ) to p * as in X v k (t 0 , x 0 ) -p * (t 0 , x 0 ) dµ 0 (x 0 ) -→ k→∞ 0. ( 23
)
Proof. See Appendix C.

Unsafe Semidefinite Programs

We will approximate programs ( 18) and ( 21) through the moment-SOS hierarchy of SDPs, as reviewed in Section 2.3. In order to apply SOS methods towards convergent approximation, we require additional assumptions:

A7 The sets X 0 , X u , X are all Archimedean BSA sets.

A8 The generator L is closed under polynomials (v(t,

x) ∈ R[t, x] =⇒ Lv(t, x) ∈ R[t, x]).
Given a degree k and a generator L obeying Assumption A8, we define the dynamics degree k as

k = ⌈ max v∈R ≤2k [t,x] deg(Lv(t, x))/2⌉. ( 24
)

Unsafe-Probability SDP

The degree-k SOS tightening of program ( 18) is

d * k (t 0 , X 0 ) = inf γ∈R γ (25a) γ -v(t 0 , x) ∈ Σ[X 0 ] ≤2k (25b) -Lv(t, x) ∈ Σ[[0, T ] × X] ≤2 k (25c) v(t, x) ∈ Σ[[0, T ] × X] ≤2k (25d) v(t, x) -1 ∈ Σ[[0, T ] × X u ] ≤2k . ( 25e
) Remark 3. Constraint (25d) absorbs the polynomial restriction of v ∈ R[t, x] ≤2k .
We require the following Lemma to ensure convergence of ( 25) to [START_REF] Ahbe | Region of attraction analysis of nonlinear stochastic systems using polynomial chaos expansion[END_REF] as k → ∞.

Lemma 5.1. Under Assumptions A1-A6, all feasible measures µ = (µ 0 , µ, µ c , µ τ ) in ( 16) are bounded.

Proof. Boundedness of nonnegative measures will be proven by the sufficient condition of compact support and finite mass. Compact support of measures in µ is ensured by A1. The initial measure µ 0 is a probability distribution (⟨1, µ 0 ⟩ = 1) by constraint (16c). The sum µ c +µ τ likewise has mass 1 by the Liouville constraint (16b), when passing in the test function v = 1. Given that µ c and µ τ are each nonnegative measures, it holds that they both have finite mass (upper-bounded by 1). Lastly, assignment of v = t to (16b) with Lt = 1 results in ⟨1, µ⟩ = ⟨t, µ c + µ τ ⟩ ≤ T . All measures have bounded masses and compact supports, and therefore are bounded.

Theorem 5.2. Under assumptions A1-A8, program (25) will converge to [START_REF] Ahbe | Region of attraction analysis of nonlinear stochastic systems using polynomial chaos expansion[END_REF]

with lim k→∞ d * k (t 0 , X 0 ) = d * (t 0 , X 0 ).
Proof. This convergence will hold by Corollary 8 of [START_REF] Tacchi | Convergence of Lasserre's hierarchy: the general case[END_REF], because all sets are Archimedean (A7), dynamics are polynomial (A8), measures are bounded (Lemma 5.1), and the objective value of ( 25) is finite (bounded below by 0).

Risk SDP

The degree-k SOS tightening for the risk-contour program in [START_REF] Chen | Safety Verification of Nonlinear Polynomial System via Occupation Measures[END_REF] is

J * k (t 0 , µ 0 ) = inf X v(t 0 , x)dµ 0 (x) (26a) -Lv(t, x) ∈ Σ[[0, T ] × X] ≤2 k (26b) v(t, x) ∈ Σ[[0, T ] × X] ≤2k (26c) v(t, x) -1 ∈ Σ[[0, T ] × X u ] ≤2k . (26d)
Remark 4. The function v(t, x) = 1 is feasible for constraints (26b)-(26d) at every degree k ∈ N. As a result, the objective in ( 26) is always upper-bounded by J * k (t 0 , µ 0 ) ≤ 1. Theorem 5.3. Under assumptions A1-A8, the risk contour program (26) will converge to [START_REF] Chen | Safety Verification of Nonlinear Polynomial System via Occupation Measures[END_REF] with [START_REF] Ahbe | Region of attraction analysis of nonlinear stochastic systems using polynomial chaos expansion[END_REF] with lim k→∞ J * k = J * . Proof. The proof of this theorem follows the same steps as in Theorem 5.2 with respect to the measure program in [START_REF] Henrion | Moment-SOS hierarchy and exit time of stochastic processes[END_REF].

Let v k (t, x) ∈ Σ[[0, T ]×X] ≤2k be the solution to [START_REF] Jasour | Risk Contours Map for Risk Bounded Motion Planning under Perception Uncertainties[END_REF] at degree-k, and let I u be the 0/1 indicator function on the unsafe set X u . Then the probability of unsafety when starting at a point x 0 ∈ X is upper-bounded by

q 1:k (x) = min(1, min k ′ ∈1..k v k ′ (t 0 , x)). (27)
Corollary 2. The sequence of functions q 1:k (x) in increasing k will converge in measure µ 0 to v(t 0 , x).

Proof. This corollary follows from Theorem 4.4, in which the sequence {v k ′ } is used to approximate v * . The minimization among all k ′ in ( 27) further sharpens the estimate of v * .

Computational Complexity

We will quantify computational complexity of the degree-k tightenings of ( 25) and ( 26) by the size of the maximal Gram matrices involved in their SOS program. In the typical case where k > k (only violated under A1-A8 when L maps every polynomial to a constant), the largest Gram matrix will occur in the Lie constraints (25c) and (26b). The Lie constraints each have n + 1 variables (t, x), so the Gram matrix size when using the monomial basis is n+1+ k k

. All other constraints have a lower degree (k rather than k), or are posed only over the n variables x. The complexity of using an interior-point method to solve the SOS programs will therefore scale based on O((n + 1) 6d ) for fixed d or O(d (n+1) ) for fixed n [START_REF] Lasserre | Moments, Positive Polynomials And Their Applications[END_REF][START_REF] Miller | Data-Driven Superstabilizing Control of Error-in-Variables Discrete-Time Linear Systems[END_REF].

Numerical Examples

MATLAB (2021a) code to reproduce all examples is available at https://github.com/jarmill/prob_ unsafe. All programs are modeled using YALMIP [START_REF] Lofberg | YALMIP : a toolbox for modeling and optimization in MATLAB[END_REF] and solved using Mosek 10 [41]. All examples will involve an initial time of t 0 = 0.

Two-Dimensional Cubic SDE

Our first demonstration analyzes safety of a cubic polynomial SDE from Example 1 of [START_REF] Prajna | Stochastic Safety Verification Using Barrier Certificates[END_REF]:

dx = x 2 -x 1 -x 2 -1 2 x 3 1 dt + 0 0.1 dW. ( 28
)
Safety of ( 28) is evaluated within the state space of X = [-2, 2] 2 until a time horizon of T = 5. The unsafe set is a moon-shaped region

X u = {x ∈ R 2 | 0.9083 2 ≤ (x 1 +0.5006) 2 +(x 2 +0.2902) 2 , 0.5 2 ≥ (x 1 -0.2) 2 +x 2
2 }. The initial set X 0 is a circle of radius R 0 and center [0.85; -0.75]. Figure 2a plots trajectories of (28) starting from X 0 (magenta region with R 0 = 0.2). Note how some sampled trajectories touch the right corner of the red half-circle X u (and are therefore unsafe). 2a plots the unsafety upper-bounding function min(1, v(0, x)) at T = 5 acquired by solving [START_REF] Sloth | Safety Analysis of Stochastic Dynamical Systems[END_REF] at order k =6 with R 0 = 0.2. The unsafe set X u is the red moon. The magenta circle is the boundary of the initial set X 0 . Note how the probability estimate is sharper in the region surrounding X 0 , as compared to the sea of Prob = 1 saturation away from X 0 .

In contrast, Figure 2b solves the order-6 SOS tightening of program [START_REF] Chen | Safety Verification of Nonlinear Polynomial System via Occupation Measures[END_REF] to produce a risk map v valid in X = [-2, 2] 2 . This risk map results in v(0, [0.85; -0.75]) = 0.4366, which is looser than the order-6 probability bound of 0.3696 from Table 1. However, Figure 2b produces an interpretable visualization of risk across the entire set X.

(a) solution of [START_REF] Sloth | Safety Analysis of Stochastic Dynamical Systems[END_REF] (b) solution of ( 26)

Figure 2: Risk levels of the unsafety upper-bound function v(0, x) at t = T = 5 and order k =6 for SDE [START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF] initialized in a disk X 0 (purple boundary) with unsafe region X u (red moon).

Adjustment of Time Horizons

In the second experiment, we consider a linear SDE with dynamics

dx = [-x 2 ; x 1 ]dt + [0; 0.1]dW (29)
The unsafe set for this example is a half-circular region with X u = {x | 0.5 2 ≥ (x + 0.5) 2 + (x 2 + 0.75) 2 , x 2 ≥ -0.75}. The plots in Figure 3 are obtained by solving (26) at order k = 6 for time horizons of T ∈ {1, 3, 5} and a state set of X = [-2, 2] 2 . As the time horizon T increases, so too does the area of X with high probabilities of unsafety. The growth of the high-probability regions occurs approximately in a clockwise pattern from X u , matching the interpretation of (29) as simple harmonic motion perturbed by diffusion in the vertical coordinate.

Discrete-time example

This subsection will focus on a discrete-time stochastic process involving a time-step of τ = 1. The law of this stochastic process with parameter λ ∈ R is

x + = -0.3x 1 + 0.8x 2 + x 1 x 2 λ/4 -0.9x 1 --0.1x 2 -0.2x 2 1 + λ/40 , (30)
in which λ is i.i.d. sampled according to a unit-normal distribution at each time (λ t ∼ N (0, 1)).

We evaluate the safety of (30) with respect to the state set X = [-1.5, 1.5] 2 , the time horizon of T = 10, and the half-circle unsafe set of [START_REF] Hilbert | Über die Darstellung definiter Formen als Summe von Formenquadraten[END_REF] evolve starting at an circular initial set with radius R 0 and center [-1; 0]. Figure 4 plots 5,000 sampled trajectories (blue dots) with respect to the magenta circular initial set and the red half-circle unsafe set. Some of the sampled trajectory points fall inside the red half-circle, corresponding to unsafety when beginning in X 0 .

X u = {x | 0.4 2 ≥ (x 1 -0.8) 2 + (x 2 -0.2) 2 , x 1 + x 2 ≥ 1}. Stochastic trajectories of
Table 2 reports probabilities of unsafety for (30) found by solving (25) at initial radii R 0 = 0 and R 0 = 0.4. To improve numerical conditioning, we normalize the time-steps from (τ, T ) = (1, 10) to (τ, T ) = (0.1, 1) without affecting the autonomous dynamics in [START_REF] Hilbert | Über die Darstellung definiter Formen als Summe von Formenquadraten[END_REF].

Figure 5a plots risk contours found by solving (25) at order 6 for R 0 = 0.4. Figure 5b plots risk contours v(0, x) acquired from (26) at degree k =6 and time t = T = 10. The returned risk map has an evaluation of v(0, [-1; 0]) = 0.4915, as compared to the (25) estimate of 7.052 × 10 -4 at R 0 = 0.

Conclusion

This paper presents a method to analyze the probability of unsafety for stochastic processes by forming LPs in occupation measures. These LPs are upper-bounded by the moment-SOS hierarchy of SDPs, yielding a convergent sequence of bounds to the true probability of unsafety. Modification of the objective leads to the development of visually interpretable risk contours for use in analysis and motion planning. The presented method can be used for any stochastic process satisfying assumptions A1-A8 (compactness, described by generators, having a polynomial description). Future work involves finding redundant constraints to refine and improve the unsafe-probability bounds. Other work includes extension to Conditional Value-at-Risk-minimizing stochastic control methods.

C Risk Contour Convergence

This appendix will provide a proof of Theorem 4.4.

We first notice that 0 ≤ p * (t 0 , x 0 ) = sup{⟨1, µ p ⟩ | µ p + µ c = δ (t0,x0) + L † µ} = inf{v(t 0 , x 0 ) | Lv ≤ 0 ; v ≥ I Xu } by applying Theorem 4.2 (using A1-A6) to the case where µ 0 = δ x0 (to have strong duality and thus access the function interpretation of the problem). We remark that our v k is feasible for the dual problem on functions [START_REF] Chen | Safety Verification of Nonlinear Polynomial System via Occupation Measures[END_REF], so by optimality of p * (t 0 , x 0 ) for each x 0 ∈ X 0 , one has v k (t 0 , x 0 ) ≥ p * (t 0 , x 0 ). Moreover, as J * (t 0 , µ 0 ) ≤ 1 < ∞ by Lemma 4.1, each v k (t 0 , •) is µ 0 -integrable, and thus the function v ⋆ (t 0 , •) is µ 0 -integrable as well. It remains to prove that their difference v k (t 0 , •) -p * (t 0 , •) ultimately vanishes in a µ 0 -sense. This vanishing difference is proven by contradiction: suppose that there exists a η > 0 such that for all k ∈ N, there is a φ(k) ≥ k with ⟨v φ(k) (t 0 , •) -p * (t 0 , •), µ 0 ⟩ ≥ η (which is the negation of ( 23)). Then, one has ⟨p * (t 0 , •), µ 0 ⟩ ≤ ⟨v φ(k) (t 0 , •), µ 0 ⟩ -η -→ k→∞ J * (t 0 , µ 0 ) -η.

However, as we already noticed, p * (t 0 , x 0 ) is feasible for [START_REF] Chen | Safety Verification of Nonlinear Polynomial System via Occupation Measures[END_REF], so this is a contradiction with optimality of the value J * (t 0 , µ 0 ) (given that J * -η < J * ). Such a contradiction proves [START_REF] Miller | Bounding the Distance to Unsafe Sets with Convex Optimization[END_REF], and therefore proves Theorem 4.4.

Theorem 4 . 3 .Theorem 4 . 4 .

 4344 Programs (20) and[START_REF] Chen | Safety Verification of Nonlinear Polynomial System via Occupation Measures[END_REF] obey strong duality m * (t 0 , µ 0 ) = J * (t 0 , µ 0 ) under A1-A6.Proof. This strong duality holds by slight modification of the method used in Theorem 3.5 and Appendix A. Let {v k } k≥1 be a sequence of solutions to (21b)-(21e) with ⟨v k

Figure 1 :

 1 Figure 1: Trajectories of SDE (28) (cyan) initialized in a disk X 0 (purple boundary), hitting an unsafe region X u (red moon).

Figure 3 :

 3 Figure 3: Risk contours for the linear SDE in (29) in increasing time horizon T w.r.t. unsafe region X u (red half-disk)

Figure 4 :

 4 Figure 4: Trajectories of discrete-time (30) (cyan) initialized in a disk X 0 (purple boundary), hitting an unsafe region X u (red half disk).

Table 2 :

 2 Unsafe probability upper-bounds for system (30) 1.00 0.1569 0.0103 1.871e-3 7.052e-4 R 0 = 0.4 1.00 1.00 0.9801 0.7054 0.5225 0.4017 (a) solution of (25) (b) solution of[START_REF] Jasour | Risk Contours Map for Risk Bounded Motion Planning under Perception Uncertainties[END_REF]

Figure 5 :

 5 Figure 5: Risk levels of unsafety upper-bound function v(0, x) at t = T = 10 and order k =6 for SDE (30) initialized in a disk X 0 (purple boundary) with unsafe region X u (red moon).

Table 1

 1 reports the probability of unsafety computed by program (25) for a circular initial set of radius R 0 .

Table 1 :

 1 Unsafe probability upper-bounds for system[START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF]

	order	1	2	3	4	5	6
	R 0 = 0	1 0.9442 0.6072 0.4805 0.4068 0.3696
	R 0 = 0.2 1 0.9736 0.7943 0.7091 0.6390 0.6132
	Figure

Acknowledgements

The authors would like to thank Roy S. Smith, Niklas Schmid, and the Automatic Control Laboratory at ETH Zürich for their advice and support.

J. Miller and M. Sznaier were partially supported by NSF grants CNS-1646121, ECCS-1808381 and CNS-2038493, AFOSR grant FA9550-19-1-0005, and ONR grant N00014-21-1-2431. J. Miller was partially supported by the Swiss National Science Foundation under NCCR Automation, grant agreement 51NF40 180545.

1

A Strong Duality

This appendix will prove strong duality between ( 16) and ( 18) by using arguments from Theorem 2.6 of [START_REF] Tacchi | Moment-sos hierarchy for large scale set approximation. application to power systems transient stability analysis[END_REF].

We collect the variables in ( 16) and ( 18) respectively into,

Variable spaces related to µ are,

with nonnegative subcones of

Corresponding spaces to ℓ are

Following the notation of [START_REF] Tacchi | Moment-sos hierarchy for large scale set approximation. application to power systems transient stability analysis[END_REF], we write that Y + = {0 Y } and Y ′ + = Y ′ . We note the containments of µ ∈ X + , ℓ ∈ Y ′ , and also that (X , X ′ ) form a pair of topological dual spaces under Assumption A1. The topologies for the spaces X and Y ′ are the weak-* topology and sup-norm-bounded weak topology respectively.

An adjoint pair of affine maps A and A ′ may be defined as

with vectors describing the cost and constraint terms as

Pairings with the vectors in (37) satisfy

⟨ℓ, b⟩ = γ.

Problem ( 16) may be expressed as

Similarly, the function LP in [START_REF] Ahbe | Region of attraction analysis of nonlinear stochastic systems using polynomial chaos expansion[END_REF] can be expressed as

Sufficient conditions to prove strong duality between ( 39) and (40) (by Theorem 2.6 of [START_REF] Tacchi | Moment-sos hierarchy for large scale set approximation. application to power systems transient stability analysis[END_REF]) are:

B Superlevel Approximation of Risk Map

This appendix provides the proof of Theorem 4.2.

B.1 Supermartingale Property

We first introduce the notion of supermartingales, and then note that any feasible v for (21b)-( 21e) is a supermartingale.

Proposition B.1. Any v that satisfies (21b) and (21e) is a supermartingale for the process L [START_REF] Rogers | Diffusions, Markov Processes, and Martingales[END_REF].

Proof. Let {X t } be the state-dependent distribution for a trajectory of the stochastic process L in (1b), respecting exit time t * ∧ τ X . This set of probability distributions satisfies the Martingale property of (7) (discrete-time) or ( 8) (continuous-time) as appropriate. Given that Lv ≤ 0 from ( 21), it holds that

The following lemma uses the previous supermartingale criterion from Proposition B.1 in combination with Doob's supermartingale inequality to provide a probability bound on v. Lemma B.2 (based on Lemma 6 of [START_REF] Prajna | Stochastic Safety Verification Using Barrier Certificates[END_REF]). Let v(t, x) ∈ C be a nonnegative function over [t 0 , T ] × X and also form a supermartingale with respect to L and X t . For a value λ ≥ 0, an initial point x 0 ∈ X, and a L-trajectory x(t) starting from x 0 at time t 0 , the following inequality holds:

Lemma B.2 modifies Lemma 6 of [START_REF] Prajna | Stochastic Safety Verification Using Barrier Certificates[END_REF] by stopping at time T and allowing for time-dependent functions v(t, x).

B.2 Superlevel Property

Theorem 4.2 can now be proven based on arguments from Theorem 7 of [START_REF] Prajna | Stochastic Safety Verification Using Barrier Certificates[END_REF]. By constraint (21d), the unsafe set X u is inside the 1-superlevel set of v(t, x) at every time t ∈ [t 0 , T ]. This superlevel relation implies that