domain: stat math sdu.ocean

text:
Victor Trappler
email: victor.trappler@univ-grenoble-alpes.fr

Élise Arnaud

Arthur Vidard

Laurent Debreu

Robust calibration of numerical models based on relative regret

Keywords: Calibration, Robust optimisation, Relative-regret, Shallow-water equations

Classical methods of calibration usually imply the minimisation of an objective function with respect to some control parameters. This function measures the error between some observations and the results obtained by a numerical model. In the presence of uncontrollable additional parameters that we model as random inputs, the objective function becomes a random variable, and notions of robustness have to be introduced for such an optimisation problem.

In this paper, we are going to present how to take into account those uncertainties by defining the relative-regret. This quantity allow us to compare the value of the objective function to its best performance achievable given a realisation of the random additional parameters. By controlling this relative-regret using a probabilistic constraint, we can then define a new family of estimators, whose robustness with respect to the random inputs can be adjusted.

Introduction

Numerical models are widely used to study or forecast natural phenomena and improve industrial processes. However, by essence models only partially represent reality and sources of uncertainties are ubiquitous (discretisation errors, missing physical processes, poorly known boundary conditions). Moreover, such uncertainties may be of different nature. [START_REF] Walker | Defining uncertainty: A conceptual basis for uncertainty management in model-based decision support[END_REF] proposes to consider two categories of uncertainties. On the one hand aleatoric uncertainties, coming from the inherent variability of a phenomenon, e.g. intrinsic randomness of some environmental variables. On the other hand, epistemic uncertainties coming from a lack of knowledge about the properties and conditions of the phenomenon underlying the behaviour of the system under study. The latter can be accounted for through the introduction of ad-hoc correcting terms in the numerical model, that need to be properly estimated. Thus, reducing the epistemic uncertainty can be done through parameters estimation approaches. This is usually done using optimal control techniques, leading to an optimisation of a well chosen cost function which is typically built as a comparison with reference observations. An application of such an approach, in the context of ocean circulation modelling, is the estimation of ocean bottom friction parameters, as done in [START_REF] Das | On the estimation of parameters of hydraulic models by assimilation of periodic tidal data[END_REF][START_REF] Das | Variational parameter estimation for a two-dimensional numerical tidal model[END_REF][START_REF] Boutet | Estimation Du Frottement Sur Le Fond Pour La Modélisation de La Marée Barotrope[END_REF]. Moreover, as such studies are often performed at a coastal and regional scale, those models are often designed with open boundary conditions, and require external forcings, such as tidal and wind forcing. Those are then a source of aleatoric uncertainties, that should be taken into account as in [START_REF] Brummelhuis | Parameter identification in tidal models with uncertain boundary conditions[END_REF].

If we overlook the aleatoric uncertainties by choosing a specific outcome, the optimal control of the parameters to be estimated can lead to localized optimisation [START_REF] Huyse | Free-form airfoil shape optimization under uncertainty using maximum expected value and second-order second-moment strategies[END_REF] and overcalibration, that is choosing a value that is optimal for the given situation. This value does not carry the optimality to other situations. In geophysics and especially in hydrological models, this overcalibration may lead to the appearance of abberations in the predictions as those uncertainties become prevalent sources of errors. In hydrology, uncertainties are the principal culprit of the existence of so called "Hydrological monsters" [START_REF] Kuczera | There are no hydrological monsters, just models and observations with large uncertainties![END_REF], that are calibrated models that perform really badly because some uncertainties in the modelling have been omitted, such as measurements errors. In flood modelling [START_REF] Wong | Neglecting model structural uncertainty underestimates upper tails of flood hazard[END_REF], the aleatoric uncertainty come from the structure of the model, and neglecting those uncertainties leads to underestimating hazard. In other domains as well, the aleatoric uncertainties can represent some manufacturing errors or environmental conditions such as wind direction and speed in wind turbine modelling [START_REF] Hu | Robust Design of Horizontal Axis Wind Turbines Using Taguchi Method[END_REF] or atmospheric conditions in aerospace vehicle design [START_REF] Brevault | Contributions to Multidisciplinary Design Optimization under Uncertainty, Application to Launch Vehicle Design[END_REF], thus they represent an important role in risk management.

It is then crucial to be able to take into account aleatoric uncertainties in optimisation problems. This consideration is called robust optimisation, or also robust design in [START_REF] Lelièvre | On the consideration of uncertainty in design: Optimization-reliability-robustness[END_REF], or optimisation under uncertainties [START_REF] Petrone | Robustness criteria in optimization under uncertainty, Evolutionary and deterministic methods for design, optimization and control (EUROGEN[END_REF][START_REF] Seshadri | A density-matching approach for optimization under uncertainty[END_REF][START_REF] Cook | Horsetail matching: A flexible approach to optimization under uncertainty[END_REF]. Furthermore, the distinction is sometimes made between stochastic optimisation and robust optimisation, depending on the knowledge available on the aleatoric uncertainties.

Let us denote k ∈ K the parameter to be estimated in order to reduce epistemic uncertainties. We assume that the aleatoric uncertainties can be modelled as a random vector U whose sample space is U. The probability measure associated with U is P U , and its density, if it exists, is p U . This distribution is assumed to be known, and that it is possible to sample from it. This choice is motivated by the fact that in various applications, the aleatoric uncertainties come from expert knowledge, empirical observations, or the knowledge acquired using other models, for instance by ensemble assimilation.

However, if it is not the case, [START_REF] Ning | Optimization under uncertainty in the era of big data and deep learning: When machine learning meets mathematical programming[END_REF] provides a comprehensive review of optimisation under uncertainty, especially on the modelling and sampling of the aleatoric uncertainties.

Since the source of the aleatoric uncertainties is considered external, the choice of k does not have any influence on the distribution of the random variable U, and therefore the aleatoric and epistemic uncertainties are assumed independent.

The cost function J(k, U) is a random variable in this context. It is most often defined as the squared norm of a given function G(k, U)

J(k, U) = 1 2 G(k, U) 2 (1)
For instance, in data assimilation, J describes a distance between the output of the numerical model and given observed data, plus generally some regularization terms. An example of such a function will be treated Section 4.

For a practical purpose, we assume that for every realisation u ∈ U of U, finding k u ∈ K that minimises the cost function k → J(k, U = u) is a well-posed problem, and that the optimum is unique. Additionally, the following assumptions are made: the cost function is strictly positive, and ∀ k ∈ K, the random variable J(k, U) has finite firstand second-order moments.

In this paper, we aim at finding k a robust estimator of k. The definition of robustness differs depending on the context in which it is used. Indeed, one definition of the robustness of an estimate is a measure of the sensibility of said estimate to outliers [START_REF] Huber | Robust statistics[END_REF]. This leads to the introduction of robust norms in data assimilation [START_REF] Rao | Robust data assimilation using $L_1$ and Huber norms[END_REF]. In a Bayesian framework, robustness may refer to the sensitivity to a wrong specification of the priors [START_REF] Berger | An overview of robust Bayesian analysis[END_REF].

As the distribution of U is assumed known, we are not going to consider distributional robustness, which deals with uncertainties on the distribution of U [START_REF] Rahimian | Distributionally Robust Optimization: A Review[END_REF].

Throughout this paper, robust has to be understood as satisfactory for a broad range of u, and/or as insensitive as possible to uncertainties encompassed in U.

The usual practice consists in neglecting the variability of U by setting it to an a priori value u b . In this case, k is set to the optimum

k u b of J(k, U = u b ).
There is no guarantee on the performance of k if the calibrated model is used for predictions, as the estimated value will compensate the error made by a possibly wrong specification of u b : this may be a case of overcalibration In a data assimilation context, this situation appears if u b does not properly represent the conditions in which the observations have been obtained.

Another strategy, that consists in minimising J over the joint space K × U, is not always possible or relevant. Indeed, joint optimisation increases the complexity of the optimisation, and the computed estimation of k has no reason to be robust in the end:

this kind of method does not take into consideration the intrinsic variability of the environmental variable. The worst-case approach [START_REF] Marzat | Worst-case global optimization of black-box functions through Kriging and relaxation[END_REF] is another popular method, and is based on the minimisation with respect to k of the maximum of the cost function for u ∈ U: min k max u J(k, U = u). This approach may yield over-conservative solutions, and does not take into account the random nature of U.

Accounting for the probabilistic nature of U leads to study the distribution of the random variable J(k, U), or the distribution of its minimisers k U . The latter is referred as the distribution of the conditional minimisers, notion that appeared notably in [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] and in [START_REF] Hennig | Entropy Search for Information-Efficient Global Optimization[END_REF] for a global optimisation purpose. Both approaches and related robust estimates are described in Section 2. Section 3 introduces a new class of estimators, by relaxing the constraint of optimality and defining regions of acceptability, similarly as [START_REF] Buhmann | Robust optimization in the presence of[END_REF] in discrete combinatorial problems, or [24,[START_REF] Snyder | Stochastic p-robust location problems[END_REF] in operations research. The rationale behind this relaxation is to be able to construct an estimate k which produces values of the cost function close enough from the minimal cost attainable given the configuration induced by u ∈ U, with high enough probability. This similarity will be measured using the relative regret. By adjusting either the relaxation or the confidence level, we can then define the RRE, the relative-regret family estimators. Illustration of the various described methods are given on a numerical example in Section 4.

Classical robust estimators

As mentioned before, robustness can be understood as satisfactory for a broad range of u, and/or as insensitive as possible to uncertainties encompassed in U. Under this definition, one may design robust estimators of k, either by using the moments of the cost function; or by exploiting the distribution of its minimisers.

Optimisation of the moments

Let us define µ(k) and σ 2 (k), the expected value and the variance of the cost variable for a given k as

µ(k) = E U [J(k, U)] = U J(k, u)p U (u) du (2)
σ 2 (k) = Var U [J(k, U)] = U (J(k, u) -µ(k)) 2 p U (u) du (3)
Minimising the expectation leads to the estimate k E defined by:

k E = arg min k∈K µ(k) (4)
In order to take into account the spread around the mean value, one can choose to minimise the variance, leading to k V :

k V = arg min k∈K σ 2 (k) (5)
A lot of different methods are readily available to solve these minimisation problems.

For instance, stochastic Sample Approximation [START_REF] Juditsky | Stochastic Approximation Approach to Stochastic Programming[END_REF][START_REF] Kim | A Guide to Sample Average Approximation[END_REF] is based on a finite and fixed set of samples {u i } i=1...N of U. The estimations at a given k are computed using standard Monte Carlo, resulting in the following optimisation problems:

kE = arg min k∈K N i=1 J(k, u i ) (6)
and

kV = arg min k∈K N i=1   J(k, u i ) - N j=1 J(k, u j )   2 (7)
If computationally affordable, one can perform these estimations on a regular grid on K × U. In case of expensive computer code, one can build a meta model to ease the minimisation, such as Gaussian processes [START_REF] Janusevskis | Simultaneous Kriging-Based Sampling for Optimization and Uncertainty Propagation[END_REF]. Even though k E is a reasonable choice, there is no guarantee that J(k E , U = u) will not reach catastrophic level for some u.

On the other hand, using k V will ensure stability of the cost function, but without any control of its performance. Ideally, one would want to have a small expectation and a small variance at the same time. Multi-objective optimisation is a proper tool to deal with these simultaneous and sometimes concurrent objectives, for exemple by computing the Pareto front of µ(k), σ 2 (k) as done in [START_REF] Baudoui | Optimisation Robuste Multiobjectifs Par Modèles de Substitution[END_REF].

As the computation of this Pareto front is usually hard and expensive, alternative strategies based on the minimisation of a scalarized version of the vector of objectives are often considered. Some are based on a weighted sum of the objectives, as presented in [START_REF] Grodzevich | Normalization and other topics in multi-objective optimization[END_REF] and in [START_REF] Marler | The weighted sum method for multi-objective optimization: New insights[END_REF], while some others are based on the minimisation of one of the objectives under constraints on the others, as performed in [START_REF] Lehman | Designing computer experiments to determine robust control variables[END_REF]. Both of these methods are based on an delicate choice of weights or of constraints before any computation. This choice relies heavily on a knowledge of the properties of the cost function.

To summarise, even though the notions of mean and variance are quite easily understood, getting a satisfactory estimator is not that straightforward. One could instead consider how often a particular value k is a minimiser of the cost function, leading to the notion of most probable estimate, as explained in the next subsection.

Most probable estimate

Let us consider the minimal cost attainable in each configuration brought by u. The resulting conditional minimum is denoted J * :

J * : u ∈ U -→ J * (u) = min k∈K J(k, U = u) (8)
Similarly, the function of conditional minimisers can then defined by:

k * : u ∈ U -→ k * (u) = k u = arg min k∈K J(k, U = u) (9)
Using this function, we can define the corresponding random variable K * as

K * = k * (U), (10)
and its associated density function p K * (k), that will be further referred as the density of minimisers. The mode of this density is called the Most Probable Estimate (MPE)

and is noted k MPE :

k MPE = arg max k∈K p K * (k) (11)
To give some intuition on this estimate, let us imagine that the distribution of minimisers is a dirac centered on k MPE . Then it would mean that this estimate is the minimiser of the cost function whatever the realisation of the uncertain variable, therefore optimal in all conditions. If the distribution p K * is heavily dominated by a single value, the MPE may be a good candidate for robust control. This is not so obvious in the case of a multimodal distribution.

In general, an analytical form of p K * is impossible to obtain, so an estimation pK * must be used, and its maximum computed to get the MPE. In the rest of the paper, the hat notation will indicate an estimation using numerical values of the underlying theoretical quantity. Once again, a set of samples {u i } i=1...N can be used to compute the set {k u i } i=1...N , from which one can approximate p K * . The resulting approximation and therefore its mode, is sensitive to the density estimation method. Main methods are KDE (Kernel Density Estimation) [START_REF] Silverman | Density Estimation: For Statistics and Data Analysis[END_REF], and EM (Expectation-Maximisation) [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF].

KDE is a non-parametric estimation technique based on the use of a kernel function

f . Assuming an isotropic kernel, the estimation has the following form:

pK * (k) = 1 N h dim K N i=1 f k -k u i h ( 12
)
where h is the bandwidth. In a multidimensional setting, one usually consider a kernel based on the product of 1D kernels, applied independently to all components:

f (k) = dim K j=1 f 1D (k (j) )
where k (j) is the j-th component of k. There is wide choice of available f 1D , and a popular choice is the Gaussian kernel

f 1D (x) = 1 √ 2π exp(-x 2 /2).
The EM algorithm can also be used to estimate the density, by minimising the statistical distance between the empirical distribution and a mixture of ν Gaussian densities.

The estimation has then the following form:

pK * (k) = ν i=1 π i φ(k; m i , Σ i ) (13)
where φ(•; m, Σ) is the probability density function of the normal distribution of mean m and covariance matrix Σ, and {π i } i=1...ν are the mixing coefficients.

In practice, despite the fact that those methods are well established, using them in a plug-in approach has some flaws. One of the basic assumption of density estimation is to assume that K * is a continuous random variable, hypothesis that may be violated.

Worse, the notion of mode is not well defined when the distribution of the minimisers is a discrete-continuous mixture. This may result in inconsistent estimations of kMPE when using different methods as illustrated in next subsection.

Numerical illustration

Before going further in the explanation of our approach, let us illustrate the nature of previously detailed estimators, kE , kV and kMPE on two analytical cost functions.

These functions are based on the Branin-Hoo's function, slightly modified to ensure strict positivity:

BH(x 1 , x 2 ) = 1 51.95 x2 - 5.1x 2 1 4π 2 + 5x 1 π -6 2 + 10 - 10 8π cos(x 1 ) -44.81 + 2 (14)
with x1 = 3x 1 -5, x2 = 3x 2
Using Eq. ( 14), we define the two cost functions on K × U = [0, 5] × [0, 5] as:

J BH : (k, u) → BH(k, u) (15)
J BHswap : (k, u) → BH(u, k) (16)
Even though the functions are quite similar, the asymmetric roles of k and u cause different behaviour in the estimations.

The random variable U is assumed to be uniformly distributed over U. Similarly, estimations of kMPE are depicted on Figure 2. The top row shows the contours of both functions as well as the set of conditional minimisers {k u i } 1≤i≤N in red, as defined in Eq. ( 9). The bottom row presents three approximations of the density of minimisers: the histogram in grey (bin size selected using Freedman-Diaconis from [START_REF] Freedman | On the histogram as a density estimator:L2 theory[END_REF]), the result of a kernel density estimation (KDE) with Gaussian kernels in red (using Scott's rule from [START_REF] Scott | On Optimal and Data-Based Histograms[END_REF] for bandwidth selection), and the estimation by a Gaussian mixture, calculated with the EM algorithm. The number of Gaussians has been fixed to 3, a guess based on the general shape of the histogram. Respective estimations of kMPE are also depicted using dashed lines.

For J BHswap , we can observe that those three methods give consistent results, as kMPE,KDE = kMPE,EM = kMPE,histogram ≈ 0.8. This is not the case for J BH : using Kernel density estimation (Gaussian), the estimation of k MPE is kMPE,KDE ≈ 1.5, while using the histogram and Gaussian mixture, kMPE,histogram = kMPE,EM = 3.8. This difference is explained by the accumulation of minimisers at this point: this challenges the assumption that K * is continuous. As the density estimation techniques traditionally assume this continuity, the EM algorithm fits this using a normal distribution with a very small variance, while the KDE considers a sum of Gaussian kernels of constant bandwidth, located at the same point. This particular problem highlights an issue with k MPE , as its estimation is possibly sensitive to the density approximation procedure.

Instead of just considering the optimal minimisers, we introduce a bit of leeway, and look for "acceptably not optimal" parameters. This slackness takes the form of a relaxation coefficient and its choice defines a new family of robust estimators, where each one of its member carries information on its robustness through this coefficient.

Relative regret-based family of estimators

Relaxing the optimality constraint

The density of minimisers has been estimated by optimising k → J(k, U = u) over K for different realisations of u. Instead of focusing on optimal values, we propose to construct an acceptable neighbourhood in terms of performances of the cost function as well. In order to do so, we are going to introduce a relaxation coefficient α ≥ 1, and given u ∈ U, we say that k ∈ K is acceptable when J(k, U = u) ≤ αJ * (u). Using the strict positivity of the objective function, we can define the relative-regret as the ratio J(k, U = u)/J * (u), and the acceptability of k will then be tied to the maximal value taken by this ratio.

In this context, for a given k, the set R α (k) ⊆ U is defined as the set of u, for which k is acceptable:

R α (k) = {u ∈ U | J(k, u) ≤ αJ * (u)} (17)
Figure 3 details the successive steps for the construction of the set R α . First, the conditional minimisers are computed, as shown on the top plots. Afterwards, for a given level α = 1.5, the set of acceptable k can be identified for each u ∈ U, as shown on the bottom left plot. Finally, the region R α (k) is the subset of U for which k is acceptable, as represented with a vertical slice on the bottom right plot.

Introducing the random nature of U, one can define Γ α (k) as the probability that k is acceptable given α:

Γ α (k) = P U [U ∈ R α (k)] = P U [J(k, U) ≤ αJ * (U)] (18)
In other words, Γ α (k) is the probability that J(k, U) is between J * (U) and αJ * (U).

Noting that without relaxation, i.e. when α is set to This linked to the definition of the distribution of the minimisers K * . For instance, if K is a discrete set, K * is a discrete random variable, and we can rewrite

1, Γ 1 is non-zero if the set {u ∈ U | J(k, U = u) = J * (u)}
Γ 1 as Γ 1 (k) = P U [J(k, U) = J * (U)] = P U [k = k * (U)] which is the probability mass function of K * .
The motivation behind this relaxation is to take into account the local behaviour of the function around the conditional minimisers. For a given set of environmental conditions u, if the function k → J(k, U = u) is flat around its minimum k * (u), then choosing k * (u) + (for a small ) will produce a value closer to the minimum than when the function has a high curvature.

0 1 2 3 4 5 k 0 1 2 3 4 5 u u k * (u) J * (u) Conditional minimiser 0 1 2 3 4 5 k 0 1 2 3 4 5 u J(k, u) < αJ * (u) J(k, u) > αJ * (u) k Relaxation of the constraint, α = 1.5 0 1 2 3 4 5 k 0 1 2 3 4 5 u Set of conditional minimisers 0 1 2 3 4 5 k 0 1 2 3 4 5 u k R α (k) Region R α (k)
We consider in this work exclusively a multiplicative relaxation, instead of an additive one in the form J(k, u) ≤ J * (u) + β for β > 0. Both formulations of the regret: the additive J -J * and the relative J/J * , take into account the value of the conditional minimum J * (u), but relative regret allow us to scale the region of acceptable points with respect to the best situation available: if J * (u) is close to 0, the region of acceptable points will grow slowly, so at a given α, we put emphasis on k * (u). On the other hand, if J * (u) is large, it means that tuning the parameter to k * (u) will not lead to very good performances of the cost function anyway. We can then put less weight on this value k * (u) in the estimation, hence the region of acceptability that grows quickly. Finally, the relative regret is a normalized quantity, in constrast to the additive regret, which shares the same unit as J. This normalization, as well as the direct relation to the relative error

J-J * J *
allows for an interpretation of the relaxation as a fraction of error.

The choice of the relaxation constant α can be made to ensure the existence of a parameter that is "acceptable" with a certain probability. For instance, given that J > 0, Γ α (k) is increasing with respect to α for any k ∈ K. We can then focus on the smallest value of α such that Γ α reaches a certain level of confidence p ∈ [0, 1]. This leads to the definition of α p

α p = inf {α ≥ 1 | ∃k p ∈ K, Γ α (k p ) ≥ p} = inf α ≥ 1 | max k∈K Γ α (k) ≥ p (19)
Rewriting the equation above, we can express α p as the solution of the following chance constrained problem

   min q s.t. max k P U J(k,U) J * (U) ≤ q ≥ p (20)
that is the smallest α, such that there exists a particular k p ∈ K for which J(k p , U) ≤ α p J * (U) with probability p. As highlighted by the formulation of Eq. ( 20), k p and α p are the result of the optimisation of the Value-at-Risk of the random variable J(k, U)/J * (U), which is a measure of risk usually applied in the financial sector (see [START_REF] Rockafellar | Deviation Measures in Risk Analysis and Optimization[END_REF]). In this spirit, when choosing k p , the relative error of the function J will be less than α p with probability p. The maximal relative regret of the function will be α p , except for the 100(1p)% least favourable cases.

The set of maximisers of Γ αp for different p: {k p for p ∈ [0; 1]} is what we are calling the relative-regret family of estimators RRE.

Figure 4 shows examples of Γ αp for J BHswap , at different levels p, and the associated estimates of the RRE. We can see that changing the level p shifts the maximiser of Γ αp .

Since k 1 is located quite far from the conditional minimisers, it arises as a compromise when the relaxation is large enough.

0 1 2 3 4 5 k 0.2 0.4 0.6 0.8 1.0 Γ α k 1 k 0.8 k 0.6
Illustration of the probability Γ αp for different p and associated k p We have then introduced the RRE, a family of estimators among which we can choose either large level of confidence p i.e. a large robustness as we are controlling the relative regret for a large fraction of possible configurations, or we can look for almost optimal performances albeit for a more reduced number of situations. Those quantities require the evaluation and optimisation of probabilities, so solving such a problem can be quite challenging. Moreover, the choice of a member of the RRE, that is the choice of a level of confidence p is also up to the modeller.

Γ α (k), α = α 1 Γ α (k), α = α 0.8 Γ α (k), α = α 0.6

Choice and computation of the relaxation coefficient

Almost-surely bounded relative-regret

Let us first consider that we want to satisfy the chance constraint of Eq. ( 20) almost surely: this is the particular case of p = 1. Given the strict positivity of J, the choice of such a level of confidence is possible only if there exists a k such that the ratio J(k, U)/J * (U) is bounded almost surely, and α 1 is the solution of the (almost-surely) chance constrained problem:

   min q s.t. min k J(k,U) J * (U) ≤ q a.s. (21)
thus, k 1 verifies

P U J(k 1 , U) J * (U) ≤ α 1 = 1 ⇐⇒ J(k 1 , U) J * (U) ≤ α 1 a.s. ( 22
)
Under the assumption that J(k, •) is continuous for all k, this is equivalent to approaching this problem using uncertain sets to model the uncertain nature of U:

J(k 1 , u) J * (u) ≤ α 1 ∀u ∈ U (23)
it follows then that

α 1 = max u∈U J(k 1 , U = u) J * (u) = min k∈K max u∈U J(k, U = u) J * (u) (24)
This can then be linked to Savage's minimax approach [START_REF] Savage | The Theory of Statistical Decision[END_REF], which consists in a worstcase approach for the additive regret J -J * , when in Eq. ( 24) we look to minimise the worst-case scenario in terms of the ratio.

Using Sample Average Approximation (SAA), based on N i.i.d. samples of U:

{u i } 1≤i≤N ,
we can reformulate the problem Eq. ( 21) as

   min q s.t. min k J(k,u i ) J * (u i ) ≤ q for 1 ≤ i ≤ N (25)
For a solution α1 of Eq. ( 25), as SAA acts as a relaxation of the initial problem Eq. ( 21), the solution found acts as lower bound on the true value α 1 : α1 ≤ α 1 .

Moreover, this estimated value α1 can be used for the estimation of the relaxation constant for a level p close to 1. Using Clopper-Pearson intervals [START_REF] Luedtke | A sample approximation approach for optimization with probabilistic constraints[END_REF][START_REF] Thulin | The cost of using exact confidence intervals for a binomial proportion[END_REF], if α1 is a solution of Eq. ( 25), it is also a feasible solution of Eq. ( 20) at a level p = η 2 1/N with probability 1η: α1 is then a probabilistic upper bound on α

( η 2 ) 1 N . α1 ∈ α ( η 2 ) 1/N , α 1 with probability 1 -η (26)
The estimated value α1 will then be optimistic in the sense that it will underestimate the true value α 1 . We can however control this approximation using the probabilistic bound lower bound of Eq. ( 26)

Choosing a level of confidence p = 1 suffers from the same pitfall as the worst-case approach, as it may return over-conservative solutions, provided that a solution exists in the first place.

Estimation of RRE

We are now going to focus on the more general case, where we can choose: to fix p ≤ 1 and deduce α p , or to set a maximal threshold α and maximise the measure of the acceptable region, or finally to find a compromise p and α p in order to keep the latter not too large.

A first approach is to set α, and then to estimate the function k → Γ α (k) as defined in Eq. ( 17). For instance, when N samples from U are available, namely

{u i } 1≤i≤N , a possible estimate of Γ α is Γα (k) = Card u i | J(k, U = u i ) ≤ αJ * (u i ) N = 1 N N i=1 1 {J(k,U=u i )≤αJ * (u i )} ( 27
)
This expression is maximised with respect to k, giving p = max k Γα (k). We can also provide a confidence interval for the true value of Γ α (k) where k = arg max Γα .

However, unless the modeller has some precise idea of desired relaxation, doing so may lead to an unsatisfactory pairing of α and p. Indeed if α is chosen too small, the resulting p = max k Γα (k) will be also small, meaning that the cost function will have non acceptable values with high probability.

Similarly, if p is fixed, the corresponding αp is computed by searching for the smallest α satisfying max k Γα (k) = p, or equivalently, by minimizing the quantile of order p of the relative-regret, which is the ratio J/J * . Once again, if αp is too large, the cost function may not be controlled enough for the contemplated application.

Looking for a compromise between p and α would be preferable. This could be achieved by studying p → α p , and particularly its slope. If this curve presents a steep increase, the multiplicative constant α p must be increased by a large amount in order to increase the probability p by a small amount. Interesting couples (p, α p ) would then be the ones located before an abrupt increase of the slope of p → α p .

Another possibility is to model this compromise by the ratio (p/α p ), as it increases with respect to p and decreases with respect to α p . The level of confidence p ratio is then defined as the maximiser of p → p/α p .

Numerical Illustration

In this Section, we will compare the different estimators introduced previously and summarised in Table 1 As stated before, we chose to model the uncertainties as a random variable uniformly distributed on U. The bounded nature of U allows us to consider members of the RRE up to a level of confidence p = 1. From now on, kMPE is estimated using KDE with Gaussian kernels.

The smallest estimated relaxation α1 and the corresponding k1 has been computed for J BH and J BHswap , using a regular grid of 1000 × 1000 points on K × U. The contour plots of those functions can be seen in the top plots of In order to choose a satisfying level of confidence p, we are going to study p → αp and p → p/α p , as described in Section 3.2.

Γα1 (k) Γ1.5 (k) k1 {(k, u) | J(k, u) = α1 J * (u)} {(k, u) | J(k, u) = 1.5J * (u)}

Histogram of the minimisers

The plot of p → αp for J BH on Figure 6 shows what seems to be a piecewise linear For J BH , the numerical values of the robust estimators can be found in Table 2. For this particular problem, the different estimates are close to each other. Practically speaking, in order to compare the effective values taken by the objective function given an estimate k we are going to consider the functions u → J( k, U = u), that we will call "profiles of k". Those profiles are well suited for the representation of the cost function for an estimate k fixed as the uncertain variable is modelled with a 1D uniform random variable.

For J BH , the curves are plotted in Figure 7. for J BHswap as no abrupt change of slope is easily discernable and the ratio presents a unique maximum for pratio = 0.766. The numerical values of the estimations k presented in Table 3 show that, contrary to J BH , the calibrated values are more spread over K. Profiles of the different estimates of J BHswap are shown in Figure 9.

In this case, kMPE , kE and kratio present a similar behaviour. They perform very well for u > 2, especially for kMPE which is very close to the minimal value; however for u < 2, they produce high values of the function. The performances of k1 are closer to the performances of kV for this function, but it performs worse than kE and kV for 4. Robust calibration of a numerical model

Calibration of a toy numerical model

We will follow the approach described in [START_REF] Kennedy | Bayesian calibration of computer models[END_REF] in order to establish the function G described in the first section in Eq. ( 1), and the resulting cost function J.

The calibration of a numerical model is usually based on the comparison between the numerical model and some observations, during a fixed time interval [0, T ] called assimilation window. The modelled physical system can be seen as a map from U to Y, the space of observations, denoted as M o : u → M o (u), where u ∈ U is an input representing some environmental conditions. The observation mentioned above is the output of the physical system during the time-window, and is denoted by

M o (u true ) ∈ Y,
where u true ∈ U is unknown.

In addition of u, the numerical model M depends on some other input k ∈ K. This additional parametrization comes usually from the successive simplifications needed to implement a numerical model of the observed physical system. k needs to be calibrated accordingly, so that the numerical model can be used to predict the behaviour of the physical system under different operating conditions.

The misfit G is defined as the difference between the numerical model and the observation. Choosing a squared norm, the cost function J defined in Eq. ( 1) is

J(k, U = u) = 1 2 G(k, u) 2 = 1 2 M(k, u) -M o (u true ) 2 (28)

The Shallow Water equations

The model to calibrate is an implementation of the 1D Shallow Water equations, described in Eq. ( 29), where h is the height of the water column, q the discharge, and z the bathymetry, while g is the usual gravitation constant. The parameter to calibrate k is the quadratic friction term, proportional to the square of the inverse of Manning-Strickler coefficient. The environmental parameter u is the amplitude of a sine wave of period 1/ω 0 . The domain of those two parameters are K = [0.0, 1.3] and U = [0.5, 0.7].

               ∂ t h + ∂ x q = 0 ∂ t q + ∂ x q 2 h + g 2 h 2 = -gh∂ x z -kq|q|h -7/3 h(0, t) = 20.0 + 3 • sin 2πt 2 + 1.5 • u • sin 2πt ω0 ∂ x q(0, t) = 0 (29)
These equations are integrated using a finite-volume scheme on a discretized domain [0, L], up to a time T . The output of the computer code is the sea surface height h, on the center of all the volumes and at all the time-steps, that will be denoted M(k, u; ω 0 = 1.0).

In this setting, the random variable U is uniformly distributed on U.

To generate the observation, we set u true = 2/3, and define M o based on the computer model M, such that M o (u true ) = M(k true , u true ; ω 0 = 0.999). ω 0 represents here the uncontrollable error between the observations and the numerical model and will now be omitted systematically in the notation.

The true value of the bottom friction k true = (k true 1 , k true 2 , . . . , k true N vol ) is not constant over the whole domain, and is defined as

k true i = 0.2 • 1 + sin 2πx i L (30)
where x i is the center of the i-th volume. The two sources of systematic errors are the one-dimensionality of K, and ω 0 . Given this setting, there exists no couple (k, u) ∈ K × U reproducing exactly the observations, thus the cost function will always be strictly positive.

Computation of the robust estimates

As the numerical model is expensive to run, it has first been evaluated on a relatively small regular grid covering K × U, and a metamodel based on Gaussian processes is constructed using those initial evaluations. In order to better capture the locus of the conditional minimisers {(k * (u), u) | u ∈ U}, we select points maximising the PEI criterion [START_REF] Ginsbourger | Bayesian Adaptive[END_REF][START_REF] Bossek | Learning Feature-Parameter Mappings for Parameter Tuning via the Profile Expected Improvement[END_REF] in order to add points to the design that improve the metamodel accuracy close to the conditional minimisers. Afterwards, a bigger regular grid is evaluated by the metamodel once the design space has been sufficiently explored, and these computations are used to calibrate the model.

The different steps of the estimation are illustrated Figure 10, where the top plot

shows the contour plot of J, with the conditional minimisers and the RRE of level p = 1. Γαp for p = 0.7, 0.8, 0.9 and 1 Looking at p → α p and p → p/α p on Figure 11, we can see that α p evolves almost linearly for p > 0.5, and that the ratio p/α p is monotonically increasing, so that the maximal ratio is found for p = 1. The different estimates take a wide range of values, as seen on Table 4, from 0.249 to almost 1.0, while k true i ∈ [0, 0.4] for 1 ≤ i ≤ N vol . It can be noted that the calibration may lead to values outside the range given by the true ones. This can be interpreted as the fact that the calibration look to compensate for errors on ω 0 . As a basis for comparison, the global minimiser of J over K × U has also been computed, and kglobal is then obtained by discarding the u value.

On

Similarly as for J BH and J BHswap , the profiles are depicted Figure 12.

We can see that the performances of k1 and kE are very similar, but k1 has better performances when u > 0.6. When comparing with the global optimiser kglobal , k1 performs better when u < 0.625, so more than half of the time. We are now going to compare how well some of those calibrated values compare in a forecast context. The mean squared forecast error averaged over the whole space and over all the timesteps gives an indication on the overall prediction quality of the prediction given this metric, and are represented on the right of Figure 13. We can see that kglobal performs slightly better than k1 , itself performing slightly better than kE . Averaging S( k) over the time steps between T and T pred , we have an indication on the quality of the forecast in the squared sense depending on the spatial coordinate x, as seen on the left of Figure 13. Relative error of the forecast error squared compared to the best attained We can see that the mean forecast error squared, on the right side, averaged over time and space is the smallest for the conditional minimisers k * (0.7) and k * (0.65), then kglobal , while k1.0 performs slightly better than kE .

It may seem surprising that some parameters calibrated without a robustness aspect perform better than k1 and kE , but their performances are largely dependent on the choice of u and their associated conditional minimisers. Good forecasts can be achieved as well as bad ones, as shown for u = 0.5, that leads to bad forecasts.

One issue with the squared forecast error is that it will penalize strongly the forecasts that present high variability with respect to the environmental conditions. A way to deal with this is to use another metric, that takes the probabilistic nature of the forecasts into account.

Continuous Ranked Probability Score

Given the random variables M pred ( k, U) and M o pred (U), representing the probabilistic forecast and the probabilistic observations, we can define the cumulative distribution functions (CDF) F pred (•, k) and F o pred (•) The Continuous ranked probability score (CRPS) measures the squared difference between the predicted CDF F pred using a calibrated value, and the CDF of the observations F o pred .

CRPS(k) = R (F pred (ξ, k) -F o pred (ξ)) 2 dξ (33)
The left plot of Figure 14 shows the CRPS averaged over time, where x denotes the spatial coordinate and the right plot shows the value of the CRPS averaged over time and space. The difference between the squared forecast error and the CRPS is apparent when comparing the general trends shared by the different calibrated parameters. According to the squared error, the sea water height of the physical region x = 4 is not well predicted, while around x = 8, the predictions are better. On the other hand, according to the CRPS, the region around x = 8 provide worse forecasts than when x = 4 and x = 7.

Given the properties of the two metrics, we can conclude that the region around x = 4

presents a lot of variability with respect to u, for both the true model and the numerical one. However, for x ≈ 8, there is a lot less variability, as the low squared error indicates, but probably a higher bias, due to the systematic errors between the truth and the numerical model.

The numerical evaluations of the CRPS for different parameters show the same order of performances observed for the squared error: the calibrated parameters that present the best performances for forecasts according to those two metrics are k * (0.7), k * (0.65) which is very similar to kglobal , and then k1 , and kE . This presupposes to know which u to choose for the conditional minimisation, thus having a strong insight on the value of the parameters in the first place. Relative augmentation of the CRPS wrt to the best attainable

Conclusion

In this paper, we dealt with a problem of robust calibration, or in other terms, the robust minimisation of an objective function. To adress this issue, we proposed a new set of robust estimators: the family of the RRE (Relative-regret estimators), and compared it in a forecast context to some other robust estimators, for the calibration of the bottom friction of a shallow water model.

Our approach is based on the comparison of the objective function to the optimal value it can attain, given a realisation of the uncertain variable. By constraining this ratio in terms of probability, the modeller can then minimize the frequency at which the relative error of the objective function will exceed a prescribed threshold. Alternatively, at a level of confidence p, the associated threshold α p bounds the relative deviation of the objective function for the proportion p of the most favourable cases. The resulting estimate can then be assessed robust, as it corresponds to values of the objective function as close as possible to the optimal values.

Except for simple analytical examples, the exact evaluation of a member of the RRE is very expensive computer-wise, not to say impossible. Some bottlenecks appear: the computations of the conditional minimisers k * (u) and the minimum J * (u), that require very local exploration of K for every u; and the computations of the probability of overshooting a given bound that depends on J * (u) (analogous to a probability of failure), task that requires an efficient exploration of the whole sample space U.

For the calibration of the shallow water model presented in Section 4, we chose to construct a response surface based on Gaussian processes, which is then used for the extensive computations. This construction is first based on an initial design, that is enriched to better capture the locus of the conditional minimisers using the PEI criterion introduced in [START_REF] Ginsbourger | Bayesian Adaptive[END_REF].

Practically speaking, the estimation of those quantities is very crude, so the computational cost may be very expensive. As perspective, specific strategies have to be created in order to carefully select each additional point that will be evaluated, and thus reduce the amount of evaluations needed of the numerical model. This selection can be based for instance on sequential design of computer experiments using Gaussian processes [START_REF] Ginsbourger | Towards Gaussian Process-based Optimization with Finite Time Horizon[END_REF].

Finally, the models upon which this calibration procedure have been applied are very simplistic. We plan in the future to apply this to the robust calibration of the bottom friction of a realistic model of the ocean.

Figure 1 :

 1 Figure 1: The left column concerns J BH , while the right one deals with J BHswap . Contours of both functions are plotted on the top, and curves of µ(k) and σ(k) are shown on the bottom (respective scales are not displayed). Estimates kE and kV are plotted with the dashed line.

Figure 2 :

 2 Figure 2: Top Left: J BH along with conditional minimisers in red. Bottom left: Estimated densities using KDE and EM algorithm and the histogram. The dashed lines indicate the MPE found using those methods. Right: Same quantities with J BHswap .

  has non-zero measure with respect to P U . It happens when the distribution of K * presents atoms.

Figure 3 :

 3 Figure 3: Principle of the relaxation of the constraint on J BH , and illustration of Rα(k). Top plots: Computation of the conditional minimisers k * (u). Top right plot: the set (k * (u), u) of conditional minimisers is represented in red. On the bottom left plot, for a relaxation α = 1.5, and u = 3, the acceptable k are in cyan, while the non-acceptable ones are in orange, while the frontier {(k * (u), u) | J * (u) = J(k, U = u)} is in yellow. On the bottom right plot, the set Rα(k) for α = 1.5 and k = 1.5 is in green.

Figure 4 :

 4 Figure 4: Illustration of the influence of different levels p on Γα p and on kp for J BHswap . When p increases, Γα p increases as well.

Figure 5 .

 5 The frontier corresponding to the couples of points (k, u) verifying {J(k, U = u) = αJ * (u)} has been drawn on top of these contour plots, for α = α1 and an arbitrary α = 1.5 < α1 to illustrate the effect of the acceptable region when the relaxation α changes. On the bottom plots, the curves k → Γα (k) for α = α1 and α = 1.5 along with the histograms of the minimisers are represented.One can notice that the relaxation allows us to avoid the issue brought by the accumulation of the minimisers of J BH at 3.8, as opposed to the MPE and its dependence on the estimation procedure of the distribution.

Figure 5 :

 5 Figure 5: Top: J BH and J BHswap contours. The thick yellow lines the boundaries of the acceptable region defined for α1 , the thick orange are for α = 1.5. The red dashed line is the estimation k1 . Bottom: Γα for α = 1.5 and α = α1 , and estimated density of the minimisers.

Figure 6 :

 6 Figure 6: Evolution of the couples (p, αp) and corresponding ratio p/αp for J BH . The dashed line indicates the level p associated with the highest ratio

Figure 7 :Figure 8 :

 78 Figure 7: Profiles of the different estimates for J BH , corresponding to the vertical cross sections of the contour. The shaded region corresponds to the interval [J * (u), α1 J * (u)]. The profiles of kMPE and kratio coincide in this case

  the middle plot are shown the conditional moments, kE and kV and the estimated density of the minimisers. At the bottom of the Figure, Γαp are represented for different levels p. As U is bounded, p = 1 is attainable, and k1 is evaluated.

Figure 10 :

 10 Figure 10: Procedure of robust calibration for the shallow water problem. Top: contours of J, conditional minimisers and {(k, U = u)|J(k, u) = α1 J * (u)}. Middle: Conditional moments and histogram and KDE of the conditional minimisers. Bottom: Γαp for different levels p

  Ratio p/αp as a function of p

Figure 11 :

 11 Figure 11: Evolution of αp for different levels p, and ratio, for the shallow water problem

4. 4 . 1 JFigure 12 :

 4112 Figure 12: Profiles for the calibration problem, for the shallow water case

  0.65) = 0.268 k * (0.7) = 0.065

Figure 13 :

 13 Figure 13: Squared forecast error for the shallow water case, depending on the calibrated parameter. The left figure shows the squared error time-averaged as a function of the spatial coordinate x, while the right barplot shows the relative change of the mean forecast error, averaged over time and space, taken w.r.t. to the best attained performance for k * (0.7)

  0.5) = 1.01 k * (0.55) = 0.737 k * (0.65) = 0.268 k * (0.7) = 0

Figure 14 :

 14 Figure 14: CRPS computed for different calibrated parameters for the shallow water case. The left figure shows the CRPS time-averaged as a function of the spatial coordinate x, and the right figure shows the relative change of the CRPS averaged over time and space, taken with respect to the best attained CRPS for k * (0.7)

  on J BH and J BHswap .

	Definition	Related quantities Interpretation
	arg min k∈K E

U [J(k, U)] k E Long run performances arg min k∈K Var U [J(k, U)] k V Steady performances arg max k∈K p K * (k) k MPE Most probable minimiser inf {α | ∃k p ∈ K, Γ α (k p ) ≥ p} (p,

k p , α p ) Acceptable values with fixed probability p p ratio = arg max p/α p (p ratio , k ratio , α ratio ) Maximal ratio of p and α p Table 1: Robust estimators, based on a cost function J

Table 2 :

 2 Estimation performed for J BH , sorted by value

	Estimate	Value
	kV	1.371
	kp , p = 1	1.557
	kE	1.587
	kMPE	1.628
	kratio , pratio = 0.654 1.637
	kp , p = 0.90	1.797

Table 3 :

 3 Estimations performed for J BHswap , sorted by value

	Estimate	Value
	kMPE	0.606
	kratio , pratio = 0.766 1.537
	kE	1.752
	kV	2.638
	k1 , p = 1	2.798

u > 2, even though its range is designed to stay within the interval [J * (u); α1 J * (u)].

We have seen how some classical robust estimators and the RRE behave on two different analytical problems. In addition to the usual levels of confidence such as 90% and 95%, one can also settle for an ad-hoc compromise, where p maximises the ratio p/α p . When the sample space of U is bounded, a conservative solution is to set p = 1.

We are now going to see how can robust minimisation is applied on the calibration of a numerical model. from the integration of the numerical model between the time T and a time T pred : the forecast window.

Given the probabilistic nature of the environmental conditions U, the forecasts will also be probabilistic. We will then compare M o pred (U) and M pred (k, U), that are the observation, and the numerical model on the forecast window for a calibrated k. Two metrics will be computed: the squared forecast error, and the Continuous Ranked Probability Score (CRPS) [START_REF] Gneiting | Probabilistic Forecasting[END_REF]. The parameters satisfying a robust criterion that are compared are the minimiser of the expectation kE , the minimiser of the variance kV and the relative regret of level p = 1, k1 as described in the previous section. We will also feature the global minimiser kglobal , and the conditional minimisers k * (u i ) where the chosen environmental variables are {u i } 1≤i≤4 = {0.5, 0.55, 0.65, 0.7}. Those values, even though they do not meet a robustness criterion, are introduced in order to have a more precise idea on the possible performances of a deterministic version of the calibration problem. The conditional minimiser k * (0.6) has been omitted as it results in a value very similar to the minimum of the expectation.

Squared forecast error

A first simple approach to measure the forecast quality is to take the squared difference between the numerical model and the observation for two samples of U. Given two environmental conditions u and u ∈ U, the former used to run the computer simulation and the latter to generate the observations, the squared forecast error for the parameter

Averaging over both u and u defines the mean squared forecast error, defined on every point of the spatial domain, and at every time-steps. This can be done using a Monte-Carlo approximation. As u and u are i.i.d., assuming that we have a set of samples {u i } 1≤i≤Nu , the squared forecast error can then be approximated:


Claudia Riccardi

Albert Meyer

Jean-Jacques Vasseur

Irene Russo Krauss

Luigi Paduano

Rosario Oliva

Luigi Petraccone

François Morvan
email: francois.morvan@umontpellier.fr

Daniela Montesarchio
email: daniela.montesarchio@unina.it

Irene Russo Krauss

Stability Is Not Everything: The Case of the Cyclisation of a Thrombin-Binding Aptamer

come   L'archive ouverte pluridisciplinaire

In the panorama of anticoagulant agents, inhibitors of thrombin -a "trypsin-like" serine protease with fundamental roles in blood clotting, converting soluble fibrinogen into insoluble fibrin 1 -are among the most reliable and widely exploited drugs against thrombosis. The 15-mer G-rich oligonucleotide Thrombin Binding Aptamer (TBA 15 or simply TBA), carrying the sequence 5'-GGTTGGTGTGGTTGG-3', is the best-characterized aptamer of thrombin. TBA has been proposed as a valuable alternative to classical thrombin inhibitors used in clinic, such as heparin, warfarin and bivalirudin, showing severe side effects and/or suffering from narrow therapeutic windows. 2 Upon folding into an antiparallel, chair-like G-quadruplex structure, TBA can tightly and selectively bind the fibrinogen-binding exosite I of human thrombin, inhibiting its key functions in the coagulation cascade. 3 Due to suboptimal dosing profiles, TBA did not progress to advanced clinical trials, being blocked after Phase I studies. 4 Since then, a large number of TBA analogues have been synthesized carrying either backbone modifications 5 or integrated into different nanosystems -including magnetic, 6 gold 7 and silica nanoparticles 8 -showing promising pharmacokinetic profiles, however not reaching advanced in vivo studies. In this scenario, as a general strategy to improve the in vivo properties of TBA, we herein propose a cyclization approach to obtain a better performing TBA analogue, involving the covalent connection, via a proper flexible linker, of the 3'-and 5'-ends of the oligonucleotide strand. Two major benefits are expected upon cyclizing TBA: on one side, the absence of the 3' and 5' termini should sensibly protect the oligonucleotide from nuclease degradation, thus significantly prolonging its in vivo half-life; on the other, the structural preorganization of the aptamer, imposed by its cyclic backbone, should favour the G-quadruplex formation, stabilizing it and thus enhancing its target affinity. This approach has been extensively adopted in the past to improve the general properties of peptides 9 and peptidomimetics, 10 as well as peptide nucleic acids (PNA) 11 and glycomimetics, 12 but only limitedly applied to oligonucleotides 13 in general, and, to the best of our knowledge, essentially unexploited thus far on aptamers. In particular, in this work we report the design, synthesis and biophysical characterization of an unprecedented cyclic TBA analogue, here named cycTBA (Fig. 1), realized as a proof of concept to validate the efficacy of cyclization approaches applied to aptamers. cycTBA -carrying a 20 atoms linker connecting the two extremities of the oligonucleotide, which in the folded G-quadruplex conformation point in the same direction, far away from both the TT loops recognized by the protein -was successfully synthesized combining standard phosphoramidite solid phase synthesis and Cu(I)-assisted azide-alkyne cyclization protocols (Fig. 1). In detail, this TBA analogue was prepared starting from the azide-functionalized solid support 1 14 and, after oligonucleotide assembly realized using standard phosphoramidite chemistry, a pentynyl group 14 was introduced at the 5'-end by reaction with phosphoramidite 2. After full deprotection and release from solid support 3, cyclization of linear oligonucleotide 4 was performed by addition of CuSO 4 and sodium ascorbate in a K + -containing solution, necessary to promote the G-quadruplex folding of the oligonucleotide and thus ensure the spatial proximity of its ends. The resulting cycTBA was purified by HPLC and then characterized by MALDI-TOF MS analysis, as well as, the regular TBA As observed previously the linear oligonucleotide 4 displayed a retention time in RP-HPLC (15.65 min) bigger than the cyclic one 5 (11.48 min). (Fig. S1 andS2).

5'

Figure 1. Synthetic scheme for the preparation of cycTBA. In the inset, a schematic representation of a cyclic TBA structure, with a generic linker connecting the 3'-and 5'-ends of the oligonucleotide sequence is reported.

Then, to get information on its conformational behaviour, the novel cyclic analogue was studied, in comparison with the unmodified TBA, by a combined approach including several techniques, i.e. UV, CD, CD-melting, differential scanning calorimetry (DSC), gel electrophoresis and SE-HPLC. In addition, its anticoagulant activity was determined by dynamic light scattering (DLS) analysis. In order to evaluate the effects of different saline conditions on the structuring ability of the aptamer, two different phosphate buffer solutions were analysed, containing a high content of K + (10 mM KH 2 PO 4 /K 2 HPO 4 , 70 mM KCl, 0.2 mM EDTA, pH 7.0, here indicated as K + buffer) or of Na + (i.e PBS: 137 mM NaCl, 2.7 mM KCl, 10 mM Na 2 HPO 4 , 1.8 mM KH 2 PO 4 , pH 7.2, here indicated as Na + buffer), respectively mimicking the intracellular and extracellular media. First of all, cycTBA and TBA samples were analysed by UV thermal difference spectra (TDS) at 2 μM concentration in both the selected K + and Na + buffers, recording spectra at low and high temperatures. The UV spectrum of cycTBA at 15 °C, in both saline solutions, revealed the characteristic absorption of TBA, with the double-hump band between 230 and 300 nm (Fig. S3 andS4). The obtained TDS profile of cycTBA showed a positive band at 273 nm and a negative band at 295 nm, similarly to unmodified TBA, providing the typical "fingerprint" of G-quadruplex structures (Fig. S3 andS4). 15 Then, the conformational behaviour in solution of cycTBA was investigated by CD experiments at different concentrations, i.e. 2, 20 and 100 µM. The CD spectra at 15 °C of cycTBAhere shown at 2 µM conc. as a representative example -gave, in both K + and Na + buffer solutions, two positive bands, with maxima at about 295 and 246 nm and a negative band with a minimum at ca. 266 nm (Fig. 2a, dashed lines). These spectral features are consistent with those obtained for the unmodified TBA analysed under the same experimental conditions (Fig. 2a), and with those reported in the literature. 16 Only a 2-3 nm shift was observed for the minimum (266 vs. 269 nm) and the maximum (246 vs. 248 nm) at low wavelengths with respect to TBA. These data confirmed that cycTBA is able to form a chair-like antiparallel G-quadruplex structure and that cyclization does not impair the G-quadruplex structuring ability of TBA nor affects its folding topology. Remarkably, in the Na + buffer the cyclic aptamer exhibited a lower structuration degree than in the K + -rich solution, as expected considering that K + can stabilize G4 structures better than Na + . 17 CD-melting experiments on cycTBA, monitoring at 295 nm, provided a nice sigmoidal behaviour with apparent T m values of ca. 69 and 60 °C, respectively in K + and Na + buffer solutions, indicating the formation of very stable G4 structures. The measured T m values were in very good agreement, within the experimental error, at all the analysed concentrations (Fig. 2, b andc). As expected, the cyclic structure allowed a marked increase of the G-quadruplex thermal stability of cycTBA, with T m = +18 °C vs. unmodified TBA (Fig. 2, S5 and S6), analysed under the same experimental conditions at 2 µM conc. (69 vs. 51 °C and 60 vs. 42 °C, respectively in K + and Na + buffer solutions). In all cases, the heating and cooling profiles were essentially superimposable, indicating that, under the experimental conditions used (scan rate of 1 °C/min), the related denaturation/renaturation processes were reversible (Fig. S5 andS6). The melting curves were also studied by van't Hoff analysis, allowing T m values consistent with those estimated from the maxima of the first derivative plots. All the thermodynamic values obtained from the CD measurements at 2 µM conc. are shown in Table 1. CD spectra were also acquired every 5 °C during the melting and cooling experiments (Fig. S7 andS8) showing a marked reduction of the CD signals at high temperatures, overall indicating that, from a qualitative point of view, cycTBA behaved as unmodified TBA in both saline phosphate buffer solutions. Notably, cycTBA was fully denaturated at 90 °C and, after the heating/cooling processes, full recovery of its spectral features was always observed (Fig. S7 andS8). Inspection of the thermodynamic parameters for cycTBA unfolding compared with TBA (Table 1) revealed that the higher thermal stability of cycTBA in both buffer solutions was mainly due to a decrease of the unfolding entropy change. This is likely the thermodynamic consequence of the structural preorganization of the aptamer, imposed by its cyclic backbone, to form the proper Gquadruplex structure. The thermodynamics of cycTBA unfolding was further investigated by differential scanning calorimetry. A noteworthy advantage of DSC measurements, with respect to spectroscopic methodologies, is that the thermodynamic parameters for the unfolding process can be obtained from experimental curves without any model assumption, whereas van't Hoff analysis of the spectroscopic melting profiles involves assuming a two-state process. The thermograms of cycTBA in both buffer solutions showed an almost symmetric shape and were fully reversible (Fig. S9). The measured enthalpy and entropy changes (Table S1) were very close to the values determined by the CD melting experiments, suggesting that the two-state assumption for the unfolding process was correct. Notably, the T m values were in very good agreement with those obtained by CD measurements, thus providing a strong evidence of the monomolecular nature of the observed transitions. Gel electrophoresis analysis was used to further characterize cycTBA in the selected buffer solutions. Different oligonucleotide concentrations, ranging from 20 to 100 μM, were explored. In Fig. S10 a representative example of a 15 % polyacrylamide gel, run under native conditions, is reported. In comparison with TBA, also cycTBA migrated as a single band on the gel, but with a slightly increased mobility, indicating that the cyclization provided the aptamer with a more compact structure, however with similar overall size in both the selected buffers (Fig. S10a).

K + buffer T m (°C) ± 1 ΔH v.H (kJ mol -1 ) ΔS v.H (kJ mol

Remarkably, no retarded bands, attributable to large aggregates or higher order G4 structures, were detected even at the highest investigated concentrations (Fig. S10b). Therefore, it can be concluded that the cycTBA behaviour was not dependent on the sample concentration and a single unimolecular G-quadruplex structure was formed in all cases, in accordance with the CD-melting data acquired at different concentrations.

In order to further investigate the conformational behaviour of the novel cyclic TBA, also a size exclusion chromatography (SE-HPLC) study was carried out. Analysis of cycTBA at 2 μM concentration showed only one peak in both buffer solutions, with retention times similar to those found for unmodified TBA (t R ca. 9 min, Fig. S11). These data are consistent with the presence of a single, monomolecular G-quadruplex structure, under the tested conditions, for both the modified and unmodified TBA, in accordance with native PAGE experiments.

As a main requirement for potential in vivo uses, TBA analogues have to show a notably increased resistance to nuclease degradation vs. unmodified TBA. Therefore, the nuclease resistance of cycTBA was studied in PBS buffer containing 80 % fetal bovine serum, i.e. a solution mimicking the extracellular environment in which the aptamer should exert its activity -in order to evaluate the protective effects of cyclization against enzymatic digestion of the oligonucleotide. Thus in parallel experiments, cycTBA and TBA were incubated in 80 % fetal bovine serum at 37 °C and, at fixed times, aliquots of these mixtures were collected and analysed both by HPLC on a Nucleogen SAX column and gel electrophoresis analysis until disappearance of the peak or the band corresponding to the pure oligonucleotide. Inspection of the chromatograms showed for the unmodified TBA the complete disappearance of the initial peak within 5 h (Fig. S12), while cycTBA was still present even after 7 days (Fig. 3). These results, also confirmed by gel electrophoresis analysis (data not shown), nicely demonstrated that cycTBA showed a markedly enhanced nuclease resistance, with t 1/2 of 100 h vs. 0.56 h for TBA, i.e. ca. 180-fold higher serum stability than TBA (Fig. 3 andS12). The high resistance to enzymatic degradation of cycTBA is clearly attributable to the presence of the linker connecting the 3'-and 5'-oligonucleotide ends, thus not accessible to exonucleases, and is highly desirable for potential in vivo studies. The inhibitory activity towards human α-thrombin was then determined by DLS experiments, following the conversion of fibrinogen into fibrin promoted by the protein, which results in a significant increase of light scattered intensity. In Fig. S13 the normalized scattered intensity is reported as a function of time for the system containing only fibrinogen and thrombin and for those containing also TBA and cycTBA, in a thrombin:aptamer ratio of 1:5. In the first case, the steep increase of intensity confirmed the quick formation of fibrin, which, on the contrary, was reduced by addition of the aptamers, as proved by the longer lag times and the lowest slope of the curves. Notably, in the presence of cycTBA, a higher intensity increase was observed than in the presence of TBA, suggesting that the former may be a less effective inhibitor than the latter. A proper comparison of their anticoagulant efficiency was performed by calculating the ratio between the initial slope of the normalized scattered intensities in the absence and in the presence of the aptamers, 18 and indeed cycTBA showed an inhibitory effectiveness that is 55% that of TBA (5±1 vs 9±1 for cycTBA and TBA respectively).

In conclusion, cycTBA exhibited an exceptional thermal stability (ΔT m of ca. +18 °C in both Na + and K + saline conditions) and notably increased nuclease resistance (ca. 180-fold higher) compared to unmodified TBA, however associated with a thrombin inhibition activityas evaluated by DLS experiments -ca. half of that found for unmodified TBA. This is a relevant result, plausibly attributable to the increased rigidity of the G-quadruplex core induced upon cyclization of the aptamer. Thus, not always higher stability and structural preorganization provide better biological results. These data seem to indicate that what could be more relevant for TBA bioactivity is the mutual interaction between the aptamer and the protein, as also suggested by other research groups. 5a Probably once thrombin recognizes the aptamer, there is a conformational rearrangement upon complexation for best interaction, which could be somehow impaired if the aptamer is too rigidified. Taken together, our results clearly demonstrated the favourable impact of cyclization on the physico-chemical features of TBA, opening new possible routes to the development of improved TBA analogues on varying the length and chemical features of the linker used to covalently bind the oligonucleotide extremities. Ongoing studies are currently devoted to obtain a fine-tuning of the linker features, so to successfully extend the cyclization approach to other biologically relevant Gquadruplex-based aptamers.

There are no conflicts to declare.

This work was supported by the Italian Association for Cancer Research (AIRC) (IG2015 n. 17037 to D. M.). F.M. is member of Inserm.

Figure 2 .

 2 Figure2. (a) overlapped CD spectra of TBA and cycTBA (continuous and dashed lines, respectively) at 15 °C and 2 µM conc. both in K + and Na + buffer solutions (red and blue lines, respectively). (b, c) Overlapped CD-melting profiles (recorded at 295 nm, with a temperature scan rate of 1 °C/min) reporting the fraction of folded cycTBA oligonucleotide as a function of the temperature at three different conc. (2, 20 and 100 µM: green, red and blue lines, respectively) in both the selected K + (b) and Na + (c) buffer solutions, in comparison to unmodified TBA (analysed at 2 µM conc., black line).

Figure 3 .

 3 Figure3. Enzymatic resistance experiments on cycTBA in 80 % fetal bovine serum (FBS) as monitored by SAX-HPLC analysis: (a) representative chromatograms of the oligonucleotide alone (yellow line) and of the reaction mixtures in FBS analysed at 1, 3 and 7 days after incubation (blue, red and green lines, respectively); (b) chromatographic peak area expressed as percentage of the remaining area normalized to the initial one for all the analysed time points (15 total time points in 7 days). Data are reported as mean values ± SD. The dark dashed line represents the 50 % remaining area and thus the half-life (t 1/2 ) in hours.

Table 1 .

 1 Apparent T

	-1 K -1 )

m values (estimated from the maxima of the first derivative plots) and thermodynamic parameters for the unfolding process of TBA and cycTBA as determined by means of CD measurements at 295 nm in both saline phosphate buffer solutions at 2 µM concentration. Errors on thermodynamic parameters are within 10 %.