Henri Mermoz

Gildas Mazo

Regularizing nested Monte Carlo Sobol' index estimators to balance the trade-off between explorations and repetitions in global sensitivity analysis of stochastic models

Keywords: nested Monte Carlo, Sobol' index, allocation, trade-off, regularization, global sensitivity analysis, stochastic model

Sobol' sensitivity index estimators for stochastic models are functions of nested Monte Carlo estimators, which are estimators built from two nested Monte Carlo loops. The outer loop explores the input space and, for each of the explorations, the inner loop repeats model runs to estimate conditional expectations. Although the optimal allocation between explorations and repetitions of one's computational budget is well-known for nested Monte Carlo estimators, it is less clear how to deal with functions of nested Monte Carlo estimators, especially when those functions have unbounded Hessian matrices, as it is the case for Sobol' index estimators. To address this problem, a regularization method is introduced to bound the mean squared error of functions of nested Monte Carlo estimators. Based on a heuristic, an allocation strategy that seeks to minimize a bias-variance trade-off is proposed. The method is applied to Sobol' index estimators for stochastic models. A practical algorithm that adapts to the level of intrinsic randomness in the models is given and illustrated on numerical experiments.

1 Introduction (Global) sensitivity analyses of stochastic models may be challenging. Indeed, stochastic models include two sources of uncertainty: the parameters' uncertainty and the intrinsic randomness of the stochastic model. The latter can be seen as a hidden additional random input, which may challenge the definition of meaningful sensitivity indices and their efficient estimation.

Many methods have been introduced to deal with stochastic models. A first approach proposed by [START_REF] Hart | Efficient computation of Sobol' Indices for stochastic models[END_REF] considers random Sobol'-Hoeffding decompositions of stochastic model outputs and defines sensitivity indices of such models as expectations of the resulting random Sobol' indices. A second approach, widely used in applications [START_REF] Courcoul | Modelling the effect of heterogeneity of shedding on the within herd coxiella burnetii spread and identification of key parameters by sensitivity analysis[END_REF], focuses on deterministic quantities of interest (QoIs) such as conditional expectations or conditional variances. By conditioning with respect to the uncertain parameters, the aim is to smooth the intrinsic randomness out and hence to deal with deterministic functions of the uncertain parameters only, so that sensitivity analysis (SA) methods for deterministic models can be applied. A third approach includes recently developed methods [START_REF] Fort | Global sensitivity analysis and Wasserstein spaces[END_REF][START_REF] Gamboa | Sensitivity analysis in general metric spaces[END_REF][START_REF] Veiga | Kernel-based ANOVA decomposition and Shapley effects -Application to global sensitivity analysis[END_REF]) that see stochastic models as deterministic functions with values in probability distribution spaces. Various sensitivity indices are defined on such spaces in order to measure the contributions of the uncertain parameters. A fourth approach rests on meta-models [START_REF] Castellan | Non-parametric adaptive estimation of order 1 Sobol' indices in stochastic models, with an application to Epidemiology[END_REF][START_REF] Étoré | Global sensitivity analysis for models described by Stochastic Differential Equations[END_REF][START_REF] Jimenez | Nonintrusive polynomial chaos expansions for sensitivity analysis in stochastic differential equations[END_REF][START_REF] Fort | Estimation of the Sobol' indices in a linear functional multidimensional model[END_REF][START_REF] Zhu | Global sensitivity analysis for stochastic simulators based on generalized lambda surrogate models[END_REF][START_REF] Janon | Asymptotic normality and efficiency of two Sobol' index estimators[END_REF][START_REF] Panin | Risk of estimators for Sobol' sensitivity indices based on metamodels[END_REF]. As a remark, let us note that that estimation of expectations of functions of conditional moments can be performed with specific methods, such as semi-parametric (da [START_REF] Da Veiga | Efficient estimation of sensitivity indices[END_REF][START_REF] Da Veiga | Efficient estimation of conditional covariance matrices for dimension reduction[END_REF] or multilevel Monte Carlo [START_REF] Mycek | Multilevel Monte Carlo covariance estimation for the computation of Sobol' Indices[END_REF] methods. Estimation in high-dimensional deterministic models is considered in de Castro and [START_REF] De | Randomized pick-freeze for sparse Sobol' indices estimation in high dimension[END_REF].

In the first three approaches described above, estimation relies on two nested Monte Carlo loops. Indeed, not only should the model be evaluated at many points in the input space (it is said that the input space is explored), but also the model should be run several times (it is said that the model is "repeated") at each of those explorations to estimate conditional expectations. In the first approach, these repetitions are performed when approximating the expectations of the random Sobol' indices. In the second and third approaches, model outputs are repeated when estimating the QoIs and the probability distributions, respectively. The larger the number of explorations and the number of repetitions, the better the sensitivity index estimates. However, in practice, models may be complex and computationally expensive. Therefore, getting a large number of runs may be impossible, and, therefore, finding a compromise between explorations and repetitions under the constraint of a computational cost or precision of estimates is important. For instance, [START_REF] Mazo | A trade-off between explorations and repetitions for estimators of two global sensitivity indices in stochastic models induced by probability measures[END_REF] proposed a choice of the number of explorations and the number of repetitions based on the minimization of some bound of the so-called mean ranking error of the estimators. This error measures the gap between the ranks of the theoretical indices and those of the estimators. However, a small mean ranking error does not necessarily imply that estimations are close to their theoretical values.

A nested Monte Carlo estimator is an estimator built from two nested Monte Carlo loops for estimating an expectation of a function of a conditional expectation. The compromise between the sizes of the outer and inner loops, coinciding, respectively, with the number of explorations and repetitions in SA, is a bias-variance trade-off. This question is sometimes refered to as the question of optimal allocation. Typically, the mean-squared error (MSE) is approximated and minimized to find the "best" asymptotic trade-off [START_REF] Rainforth | On Nesting Monte Carlo Estimators[END_REF][START_REF] Giles | Multilevel nested simulation for efficient risk estimation[END_REF][START_REF] Giorgi | Weak Error for Nested Multilevel Monte Carlo[END_REF][START_REF] Hong | Estimating the mean of a non-linear function of conditional expectation[END_REF]. More precisely, if T denotes the number of model runs available to the practitioner, and n = T 1-η and m = T η , η ∈ [0, 1), denote the sizes of the outer and inner loops, respectively, then it is desired to know which value of η would lead to the optimal convergence rate of the MSE.

Sobol' index estimators are functions of nested Monte Carlo estimators. It is, therefore, natural to ask whether the MSE of functions of nested Monte Carlo estimators can be bounded, too. The answer will depend on the regularity of the function in question. For instance, if the function is twice continuously differentiable with a bounded Hessian matrix, then standard Taylor expansion techniques can be used. However, if the Hessian is unbounded, as it is the case of Sobol' index estimators, then the question is more challenging and remains open, to the best of our knowledge.

We note that the question of optimal allocation of available computational resources is shared more globally by research fields interested in stochastic simulators (Chen andZhou, 2014, 2017;[START_REF] Binois | Practical Heteroscedastic Gaussian Process Modeling for Large Simulation Experiments[END_REF][START_REF] Binois | Replication or Exploration? Sequential Design for Stochastic Simulation Experiments[END_REF].

In this paper, we introduce a regularized version of the nested Monte Carlo Sobol' index estimator and bound its MSE. Based on a heuristic, the bound is simplified and a minimization is carried out to balance the numbers of explorations and repetitions. A practical algorithm that adapts to the model intrinsic randomness is proposed and illustrated on numerical experiments.

This paper is organized as follows. Sobol' indices for stochastic models and their estimation based on nested Monte Carlo methods are presented in Section 2. The regularization method, as well as the challenges of bounding the MSE of functions of nested Monte Carlo estimators, are introduced in Section 3. Section 4 is an application to Sobol' index estimators. Practical algorithms are presented and illustrated in Section 5. A conclusion closes the paper.

Sensitivity analysis of stochastic models

A stochastic model with inputs X = (X 1 , . . . , X p ) ∈ R p and output Y ∈ R is modeled as a function f of X and some collection of random variables, denoted by Z, independent of X such that

Y = f (X, Z) . (1)
The stochasticity of the model originates from Z since the output of the model evaluated at an input X = x is a random variable f (x, Z). The distribution of Z is generally unknown.

In the context of SA, a way to deal with stochastic models consists in carrying out SA for deterministic models given by deterministic QoIs. This allows to switch from a stochastic model to some deterministic models for which many SA methods are studied in the literature.

We consider QoIs of the form

Q(X) = E [ϕ(X, Z) | X] , (2)
where ϕ(X, Z) is a function of X and Z. For instance, if

ϕ = f then Q(X) is the conditional expectation of the model and if ϕ(X, Z) = (f (X, Z) -E [f (X, Z) | X]) 2 then Q(X) is the conditional variance, two common choices in practice.
If u is a subset of {1, . . . , p}, denote by X u the group of inputs {X i , i ∈ u} and X ∼u the group of inputs {X i , i ∈ u}. The Sobol' and total indices of the input vector X u with respect to the function Q are defined as

S u = Var (E [Q(X) | X u ]) Var (Q(X)) (3)
T u = 1 - Var (E [Q(X) | X ∼u ]) Var (Q(X)) = 1 -S ∼u . (4)
The Sobol' index S u (and hence T u ) can be expressed in terms of a function g linking the components of some parameter vector. Let X be an independent copy of X, independent of Z. Denote by X ∼u the subvector of X whose components are those of X not indexed by u. (For instance, if p = 4 and u = {1, 4} then

X ∼u = ( X 2 , X 3 ).) If θ = (θ 1 , θ 2 , θ 3 ) with θ 1 = E(Q(X) 2 ), θ 2 = E(Q(X)) and θ 3 = E(E [Q(X) | X u ] 2 ) = E(Q(X)Q( X ∼u , X u )) = E(Q(X ∼u , X u )Q( X ∼u , X u )) then S u = g(θ) = θ 3 -θ 2 2 θ 1 -θ 2 2 .
The construction of an estimator of S u boils down to the construction of estimators of the three quantities θ 1 , θ 2 and θ 3 . To do this, nested Monte Carlo sampling is widelyused. Let {X (i) ; i = 1, . . . , n} and { X (i) ; i = 1, . . . , n} be independent Monte Carlo samples from the law of X. For each i, denote by X (i) u the subvector of X (i) whose components are those of X (i) indexed by u. Likewise, denote by X (i) ∼u the subvector of X (i) whose components are those of X (i) not indexed by u, and denote by X (i) ∼u the subvector of X (i) whose components are those of X (i) not indexed by u. An estimator of S u is given by

S u = g( θ) = θ 3 -θ 2 2 θ 1 -θ 2 2 (5)
where

θ 1 = 1 n n i=1 Q m (X (i) ) 2 θ 2 = 1 n n i=1 Q m (X (i) ) θ 3 = 1 n n i=1 Q m (X (i) ) Q m ( X (i) ∼u , X (i) u )      (6)
and

Q m (X (i) ) = 1 m m k=1 ϕ X (i) , Z (i,k) Q m ( X (i) ∼u , X (i) u ) = 1 m m k=1 ϕ( X (i) ∼u , X (i) u , Z (i,k) );
here the objects {Z (i,k) , Z (i,k) ; k = 1, • • • , m; i = 1, • • • , n}, are independent and identically distributed random variables, independent of {X (i) , X (i) ; i = 1, . . . , n}, representing the randomness of the user's model. If m is fixed and n → ∞, then

√ n S u -S u 1 - E Var(ϕ(X, Z)|X) E Var(ϕ(X, Z)|X) + m Var E(ϕ(X, Z)|X)
converges to a centered normal distribution with some variance σ 2 m depending on m.

If, moreover, m → ∞ such that √ n/m → 0 then √ n( S u -S u )
goes to a centered normal distribution with variance lim m→∞ σ 2 m . For more details, see [START_REF] Mazo | A trade-off between explorations and repetitions for estimators of two global sensitivity indices in stochastic models induced by probability measures[END_REF]. A theoretically-guided choice for m and n that penalizes bad rankings of the sensitivity index estimates S 1 , . . . , S p was given in [START_REF] Mazo | A trade-off between explorations and repetitions for estimators of two global sensitivity indices in stochastic models induced by probability measures[END_REF].

Nested Monte Carlo estimation

We say that an estimator θ is a nested Monte-Carlo estimator if

θ = 1 n n i=1 φ 1 1 m m k=1 ϕ 1 (ξ (i) , ζ (i,k) ) • • • φ p 1 m m k=1 ϕ p (ξ (i) , ζ (i,k) ) , (7)
where φ 1 , ϕ 1 , . . . , φ p , ϕ p are measurable functions and the sets of random vectors {ξ (i) , ζ (i,k) , k = 1, . . . , m}, i = 1, . . . , n, are independent and identically distributed. The estimator ( 7) is called a nested Monte Carlo estimator because it is built from two nested Monte Carlo loops: an outer loop is used to simulate the random vectors ξ (1) , . . . , ξ (n) , and, for each ξ (i) , an inner loop is used to simulate ζ (i,k) , k = 1, . . . , m. The random vectors {ζ (i,k) } and {ξ (i) } are assumed to be independent. With this definition, each of the three estimators θ 1 , θ 2 , θ 3 in ( 6) is a nested Monte Carlo estimator: the outer loop simulates the "explorations" X (i) , X (i) , i = 1, . . . , n, and the inner loop simulates the "repetitions" ϕ(X (i) , Z (i,k) ), ϕ(

X (i) ∼u , X (i)
u , Z (i,k) ), k = 1, . . . , m. If, for each i = 1, . . . , n and conditionally on ξ (i) , the random variables ζ (i,k) are independent and identically distributed, and if

θ := E φ 1 E ϕ 1 (ξ (1) , ζ (1,1) ) | ξ (1) • • • φ p E ϕ p (ξ (1) , ζ (1,1) ) | ξ (1)
, (8) then we have

E θ -θ 2 = 1 n Trace (Σ m ) + b m 2 , ( 9
)
where Σ m is the variance-covariance matrix of

φ 1 1 m m k=1 ϕ 1 (ξ (1) , ζ (1,k) ) • • • φ p 1 m m k=1 ϕ p (ξ (1) , ζ (1,k) ) .
The MSE in Equation ( 9) is a sum of two terms: a variance term which involves the ratio 1/n and a bias term that is function of m only.

Assumption 1. b m → 0 and Σ m → Σ as m → +∞ Assumption 1 is fulfilled by Nested Monte Carlo estimators provided that the functions φ 1 , . . . , φ p in (7) have good properties such as boundedness, Lipschitz continuity or boundedness of derivatives [START_REF] Giorgi | Limit theorems for weighted and regular Multilevel estimators[END_REF][START_REF] Giorgi | Weak Error for Nested Multilevel Monte Carlo[END_REF]. In the context of global sensitivity analysis for stochastic model, [START_REF] Mazo | A trade-off between explorations and repetitions for estimators of two global sensitivity indices in stochastic models induced by probability measures[END_REF] showed that Assumption 1 is satisfied by θ 1 , θ 2 , θ 3 (see Equation ( 6)) associated with Sobol' index estimators.

If Assumption 1 is satisfied, it holds that lim n,m→+∞ E θθ 2 = 0 and thereby θ converges in quadratic mean to θ. If, moreover, the rate at which b m goes to zero is known, then optimal rates for n = T 1-η and m = T η , η ∈ [0, 1), can be calculated in terms of the total number of runs T . For instance, in the case p = 1 and φ 1 is a smooth function with uniformly bounded third derivative, [START_REF] Hong | Estimating the mean of a non-linear function of conditional expectation[END_REF] showed that

b m = O(1/m) and thus E θ -θ 2 = O(1/n + 1/m 2 ), leading to O(T 2/3 ) reached for n = O(T 2/3 ) and m = O(T 1/3 ). If φ 1 is simply Lipschitz continuous, Rainforth et al. (2018) rather found that E θ -θ 2 = O(1/n + 1/m), leading to a MSE of order O(T 1/2 ), ensured by the choice of n = O(T 1/2 ) and m = O(T 1/2 ).

Functions of nested Monte Carlo estimators

Let g : D → R be a continuous, non-constant function and let θ be as in ( 6). Note that the estimator g( θ) converges in probability to g(θ) but it does not necessarily hold that its MSE E (g( θ)g(θ)) 2 converges to 0. In fact, the MSE may even be infinite, depending on the properties of g and the law of θ. Since the law of θ is unknown, it is important to find verifiable conditions on g under which the MSE is kept under control for all possible laws of θ. For instance, if g is Lipschitz continuous, then there exists an constant L such that

E g( θ) -g(θ) 2 ≤ L 2 E θ -θ 2 .
This bound could be minimized as in [START_REF] Hong | Estimating the mean of a non-linear function of conditional expectation[END_REF], see above. If g is assumed to be a twice continuously differentiable such that its Hessian matrix denoted ∇ 2 g is uniformly bounded, then, up to existence of some moments of θ, the combination of a Taylor-Lagrange expansion and convexity inequalities would yield

E g( θ) -g(θ) 2 ≤ 4E ∇g(µ m ) ⊤ θ -µ m 2 + 2 (g(µ m ) -g(θ)) 2 + L ′2 E θ -µ m 4 , ( 10
)
with L ′ > 0 such that sup x∈D ∇ 2 g(x) F ≤ L ′ and • F denotes the Frobenius norm. From this upper bound, we could derive another bound by noticing that

E θ -µ m 4 = O(1/n 2 ).
This follows from the Marcinkiewicz-Zygmund inequality [START_REF] Marcinkiewicz | Sur les fonctions indépendantes[END_REF], provided that suitable moments exist. It would then hold that:

E g( θ) -g(θ) 2 ≤ 4E ∇g(µ m ) ⊤ θ -µ m 2 + 2 (g(µ m ) -g(θ)) 2 + O(1/n 2 ). (11)
Equation ( 11) highlights three terms in the right-hand side: the first term stands for a variance term, the second term represents a squared bias term and the third term is a negligible approximation error. However, in practice, in many interesting applications, the function g does not have good enough properties such as above. For instance, uniform boundedness of derivatives of g is a very strong condition in general. In order to weaken such a condition, consider the following one:

E sup λ∈[0,1] ∇ 2 g λ θ + (1 -λ)µ m 4 F = O(1) as n, m → ∞. ( 12
)
Under the condition (12), the decomposition (11) holds but asymptotically when both n and m go to infinity. However, whether or not this condition is true depends on the law of θ, and hence hard to check in practice. For instance, in the case of Sobol' index estimators defined in Section 2, this comes down to require that E[( θ 1θ 2

2 ) -α ] with α > 0 exists, whereas the probability distribution of θ 1θ 2 2 is unknown. In the next section, we propose a condition, weaker than the one in Equation ( 12), that does not put constraints on the distribution of θ.

A regularization method to address functions with unbounded Hessian matrix

The idea is to introduce a "slight perturbation" g h (with h ∈ (0, 1)) of the function g that approaches the true g as h → 0 and such that ( 12) is satisfied with g h in place of g. Introducing g h can be thought as a way to "transport" the original estimator θ to regions of D where control of moments of g h ( θ) is possible without additional conditions on the law of θ. The goal is to "get away" from certain regions of the parameter space where the Hessian of g may explode. For this, consider a family of functions {g h : D → R, h ∈ (0, 1)} such that for all h, g h is twice continuously differentiable and for all x ∈ D, lim h→0 g h (x) = g(x). Moreover, assume g h satisfies the following assumption:

Assumption 2. There exists a constant C independent of h such that, for all h ∈ (0, 1):

lim sup n,m→∞ E sup λ∈[0,1] ∇ 2 g h (λ θ + (1 -λ)µ m ) 4 F ≤ C.
The advantage of having such a family of functions is that E(g h ( θ)g h (θ)) 2 , the "perturbed MSE", could be bounded with an approximate upper bound in the form of Equation ( 11) with g = g h . A direct implication of this is that for h fixed,

lim n,m→∞ E g h ( θ) -g h (θ) 2 = 0. Further, if the regularized estimator g h ( θ) is
used for estimation of g(θ) instead of g( θ), the corresponding mean squared error is bounded in Proposition 1.

Proposition 1. Under Assumptions 1 and 2, for every h ∈ (0, 1):

E g h ( θ) -g(θ) 2 ≤ 4 (1 + p n,m (h)) V n,m (h) + B m (h) 2 , ( 13
)
where lim sup n,m→∞ p n,m (h) = 0 and V n,m (h) := E(∇g h (µ m ) ⊤ ( θ -µ m )) 2 and B m (h) := g h (µ m ) -g(θ).
Remark that B m (h) = (g h (µ m )g h (θ)) + (g h (θ)g(θ)) accounts for two types of bias. The first term in the sum stands for the bias of the estimator θ; it vanishes as soon as µ m → θ. The second term is a regularization error that stems from the use of g h instead of g. This term does not vanish asymptotically in general. However, it goes to zero as h → 0, yielding Corollary 1.

Corollary 1. We have lim h→0 lim n,m→∞ E g h ( θ) -g(θ) 2 = 0.
Note that the limits in Corollary 1 cannot be interchanged. Additional conditions on the gap E(g( θ)g h ( θ)) 2 would be needed to get the convergence of the MSE of g( θ).

Application to Sobol' index estimators

Let θ = ( θ 1 , θ 2 , θ 3 ) be as in (6). Recall that θ 1 = E(Q(X) 2 ), θ 2 = E(Q(X)), θ 3 = E(E [Q(X) | X u ] 2 ) = E(Q(X)Q( X ∼u , X u )), so that µ m = (µ m1 , µ m2 , µ m3 ) = (θ 1 +b m1 , θ 2 +b m2 , θ 3 +b m3 ) and b m = (b m1 , b m2 , b m3 ) = (EVar [ϕ(X, Z) | X] /m, 0, 0). Recall that S u = g( θ)
, where the function

g : (x 1 , x 2 , x 3 ) → (x 3 -x 2 2 )/(x 1 -x 2 2 ) (14)
is a twice-continuously differentiable function over its definition domain. But unfortunately, g is unbounded. The form of such a function makes the study of the MSE of S u difficult unless strong conditions are imposed on the output distribution of the stochastic model. This could explain why, in spite of our efforts, it is difficult to find in the literature material on the convergence in quadratic mean of such estimators. The approach introduced in Section 3 allows to bypass the unboundedness issue. We shall apply the results of Section 3 to bound the MSE of the regularized version of the Sobol' index estimator and propose an allocation method.

Throughout this section, it is assumed that E Q(X) 16 < +∞. We consider the family of function g h (x) := g(x + hu) with u = (1, 0, 0). In order to provide an upper bound for E(g h ( θ)g(θ)) 2 as in Proposition 1, it is necessary to fulfill Assumption 1 and 2. Assumption 1 is trivially satisfied.

Theorem 1. With θ, g h and g as above, Assumption 2 holds. There is a constant C

independent of h such that, for every h ∈ (0, 1),

E g h ( θ) -g(θ) 2 ≤ C 1 n + E (Var [ϕ(X, Z) | X]) m 2 + h 2 , as n, m → ∞.
Corollary 2. With θ, g h and g as above, we have lim h→0 lim n,m→∞ E(g h ( θ)-g(θ)) 2 = 0.

As n, m → +∞, it occurs that the smaller h is, the smaller the MSE gets. Therefore, in practice, parameter h should be taken in order to approximate the true g(θ).

To balance the trade-off between the number of repetitions and explorations, we proceed with a heuristic. For h small enough, one could assume that h 2 is negligible with respect to the quantity BVT := 1/n + (E (Var [ϕ(X, Z) | X]) /m) 2 . This is, therefore, the quantity we shall be interested in for optimizing the numbers of repetitions and explorations. Let T = mn be the number of available samples to compute the Sobol' index estimator. From nested Monte Carlo theory (see beginning of Section 3), we already know that m and n must be of order T 1/3 and T 2/3 respectively so as to obtain optimal convergence rate.

However, an asymptotic order is not a specific value. We expect that the greater the variance of intrinsic randomness, the larger the value of optimal m should be. Thus, in this section, the goal consists in proposing a value of m that adapts to the intrinsic randomness of the stochastic model.

Under the constraint nm = T , the optimal convergence rate of the BVT is obtained when n opt is of order T 2/3 and m opt is of order T 1/3 . Let m opt = κT 1/3 where κ > 0. Then, n opt = κ -1 T 2/3 . Thus:

BVT = κT -2/3 + (EVar [ϕ(X, Z) | X]) 2 κ 2 T -2/3 . Coefficient κ can be chosen such that κT -2/3 + (EVar[ϕ(X,Z)|X]) 2 κ 2
T -2/3 is the smallest over κ > 0. The minimum of such a quantity is reached at κ opt = 2 (EVar [ϕ(X, Z) | X]) 2 1/3 . Therefore:

m opt (T ) = 2 (EVar [ϕ(X, Z) | X]) 2 1/3 T 1/3 . ( 15
)
Therefore, the number of repetitions suggested above ensures first that the BVT converges at optimal rate and second adapts to the model's intrinsic randomness because that m opt (T ) depends on EVar [ϕ(X, Z) | X], which quantifies the part of the total variance Var (ϕ(X, Z)) that is not attributed to the inputs X; and so, that measures the influence of the intrinsic noise of the stochastic model ϕ(X, Z). Finally, it also appears that m opt (T ) depends on both T and the function ϕ. The dependence with respect to T guarantees that m opt (T ) grows as T gets large. Besides, the dependence with respect to ϕ means that even if m opt (T ) remains proportional to T 1/3 , it also varies with respect to the chosen QoI of the stochastic model f .

5 An algorithm of estimation of Sobol' indices for stochastic models

Algorithms

This section is devoted to the practical implementation of the bias-variance trade-off strategy when performing SA for some QoI of a stochastic model. Recall that f is a stochastic model as in (1) and we are interested in carrying out SA of a QoI under the form (2), that is,

Q(X) = E [ϕ (X, Z) | X]
in order to measure the impact of l groups of inputs u i ⊂ {1, . . . , p}, i = 1, . . . , l. In other words, we are interested in estimating S u1 , . . . , S u l . We shall use at most T × (l + 1) evaluations of ϕ(X, Z).

Under the constraint nm = T , the number of repetitions m opt found in (15) depends on ρ := EVar [ϕ(X, Z) | X]. However, in practice, ρ is often unknown. So, before sensitivity index estimation, ρ needs to be estimated. Notice that ρ can be rewritten ρ = E(ϕ(X, Z)ϕ(X, Z)) 2 , which suggests the following estimation procedure. Consider r 0 i.i.d. samples of X, denoted by r0) , and generate two outputs at each sample X (i) : ϕ(X (1) , Z (1,1) ), ϕ(X (1) , Z (1,2) ) , • • • , ϕ(X (r0) , Z (r0,1) ), ϕ(X (r0) , Z (r0,2) ) . Thus:

X (1) , • • • , X (
ρ = 1 r 0 r0 i=1 ϕ(X (i) , Z (i,1) ) -ϕ(X (i) , Z (i,2) )
2 is a consistent and unbiased estimator of ρ. It appears that the estimation of ρ requires 2r 0 evaluations of the model ϕ(X, Z). However, the maximal number of evaluations is T × (l + 1). So, for index estimation procedure, at most T × (l + 1) -2r 0 model evaluations are allowed. Therefore, our strategy consists in leveraging the model outputs used to estimate ρ and then plugging and completing those outputs in order to compute sensitivity index estimates. This strategy relies on two algorithms: Algorithm 1 and Algorithm 2. Algorithm 1 enables to generate complementary outputs in addition to outputs already available after estimation of ρ. This allows to satisfy the constraint of the maximal number of model evaluations T × (l + 1). This part helps to optimize the whole estimation procedure by using the model outputs already generated. Regarding Algorithm 2, it effectively estimates indices in three steps based on pick-freeze procedure. First, it estimates ρ and thereby compute m opt and n opt = T /m opt . Then, in the second step, by relying on Algorithm 1, complementary outputs required for estimations are generated. In the final step, sensitivity index estimates are computed with respect to inputs or groups of inputs specified by the user.

Algorithm 2 Estimation of Sobol' indices

Inputs: h, r 0 , T, ϕ, w = {u 1 , • • • , u l }, X (1) , • • • , X (T ) , X (1) , • • • , X (T ) 1: for i = 1, • • • , r 0 do 2: for k = 1, 2 do 3:
Compute ϕ X (i) , Z (i,k) 4: end for 5: end for 6: Compute ρ ← 1 r0 r0 s=1 ϕ X (i) , Z (i,1)ϕ X (i) , Z (i,2) 2 7: Compute m opt according to Equation (15) 8: n opt ← T / m opt 9: Run Algorithm 1 with m = m opt , n = n opt to complete samples ϕ X (1) , Z (1,1) , • • • , ϕ X (r0) , Z (r0,2) . 10: for j = 1, • • • , l do 11:

for i = 1, • • • , n opt do 12: for k = 1, • • • , m opt do 13: Compute ϕ ( X (i) ∼uj , X (i) uj ), Z (i,k) 14:
end for 15:

end for 16: end for 17: Compute sensitivity index estimates or their regularized estimates Algorithm 2 requires: h, r 0 , T, ϕ, w and input samples. The transformation ϕ of the stochastic model is supplied as well as w the set of inputs or groups of inputs whose indices are estimated. In practice, r 0 , T and h must be chosen. We recommend to take r 0 small with respect to T so as not to waste too much of the budget in the first stage of Algorithm 2. Indeed, the estimator ρ is unbiased and consistent and often provides good estimates even for small r 0 . Regarding T , it follows T should be taken as large as possible depending on the computational cost of a run of both ϕ and the original model f . Furthermore, to ensure that the MSE has a precision ε ∈ (0, 1) with h 2 ≪ ε, T must be roughly chosen larger than ε -3/2 since the MSE is O T -2/3 . This provides approximations for practical choice of T .

Illustrations

This section presents the performance of the estimators of first-order and total indices computed by Algorithm 2 in the case of two toy stochastic models for which analytical values of indices are known: a linear model f [START_REF] Ishigami | An importance quantification technique in uncertainty analysis for computer models[END_REF]). For each value of T = nm, the estimators of Algorithm 2 are compared with two other arbitrary choices, namely, (n, m) = (T /5, 5) and (n, m) = (T 1/2 , T 1/2 ).

(X 1 , X 2 , Z) = 1 + X 1 + 2X 2 + σZ with σ > 0 and a stochastic version of the Ishigami function f ′ (X 1 , X 2 , X 3 , Z) = sin X 1 + a sin 2 X 2 + bX 4 3 (sin X 1 )Z 2 with a, b > 0 (Ishigami
These two choices above represent two different situations. The choice (n, m) = (T /5, 5) presents a case where the number of repetitions is constant and independent of T . This illustrates the situation where the bias does not get reduced so that it disturbs estimations no matter how large T is. Regarding (n, m) = (T 1/2 , T 1/2 ), it shows that the case where the variance is not sufficiently reduced.

For illustrations, the product T = mn is chosen in the set T ∈ {10 3 , 10 4 , • • • , 10 7 }.

For each choice of the couple (n, m), N = 100 replications of estimations are carried out. The tuning parameter r 0 is set to 10. Thus, 2r 0 = 20 model evaluations are used to get the estimates ρ, n opt , m opt in the first part of Algorithm 2, and then T -2r 0 ∈ {10 3 -20, 10 4 -20, • • • , 10 7 -20} model evaluations are used to get the sensitivity index estimators with n = n opt and m = m opt . For both toy stochastic models, the QoI considered is the conditional expectation so that ϕ = f or ϕ = f ′ depending on the model. For each value of T , the boxplots of sensitivity index estimates are plotted for each of the three choices.

Linear model

Let f (X 1 , X 2 , Z) = 1 + X 1 + 2X 2 + σZ where σ > 0 and X 1 , X 2 and Z are i.i.d. under the standard normal distribution. Such model includes two uncertain parameters X 1 and X 2 with respective first-order Sobol' indices S 1 = 1/5 and S 2 = 4/5. Two values of σ are considered: σ = 1 and σ = 5.

Figure 1 shows that the estimations obtained with Algorithm 2 are more efficient as T gets large because both bias and variance are efficiently reduced. Boxplots highlight that the strategies m = 5 and m = T 1/2 suffer respectively from bias and variance. Notice that in the case of the linear model under study, E(Var(f (X 1 , X 2 , Z) | (X 1 , X 2 ))) = 2σ 2 ; so the bias depends on σ. This explains why in the case σ = 5 (see Figure 2), even for large value of T , estimations resulting of the choice m = 5 do not converge to their true targets. and m = T 1/2 (in blue). The red dashed lines showed the true sensitivity index values, in particular, in this setting S 1 = 0.2 and S 2 = 0.8 -π, π]) and Z ∼ N (0, 1). The model f ′ is a modified version of benchmark function known as the Ishigami function in SA. For this model, first-order Sobol' indices of inputs X 1 , X 2 , and X 3 for the QoI E(f ′ (X 1 , X 2 , X 3 , Z) | X 1 , X 2 , X3) are respectively given by Figures 3 and 4 also reveal that estimations obtained by using Algorithm 2 are more efficient for large T . Besides, remark that the term bX 4 3 sin X 1 that multiplies the intrinsic noise term Z 2 includes b so that |bX 4 3 sin X 1 | ≤ bπ 4 . Then, b allows to control the variance of the intrinsic noise term of the model. In both cases, the strategy implemented in Algorithm 2 has better results.

T=1e+04

A stochastic Ishigami function Let f ′ (X 1 , X 2 , X 3 , Z) = sin X 1 +a sin 2 X 2 +bX 4 3 (sin X 1 )Z 2 such that with a, b > 0, X 1 , X 2 , X 3 and Z are independent with X 1 , X 2 , X 3 distributed under U ([
S 1 = 1 2 1+ bπ 4 5 2 a 2 8 + bπ 4 5 + b 2 π 8 18 + 1 2 , S 2 =
Overall, it appears that the strategy of Algorithm 2 provides better estimations and estimators it generates converge faster. In the particular case of m = 5, it is noticeable that errors do not decrease when T gets larger but rather they are quite constant. This is explained by the fact that the bias is constant since m is constant. This illustrates the importance of varying the number of repetition when the total computational budget grows. Regarding the case m = T 1/2 , it turns out that estimators do not converge with optimal rate compared to the case m = m opt due the variance part of the meansquared error. Indeed, with m = T 1/2 , the variance converges to 0 at rate T 1/2 while the squared-bias converges at rate T . Then, the global convergence rate of the MSEs is T 1/2 that is slower than the rate T 2/3 of estimators built by Algorithm 2. These two cases clearly illustrate the bias-variance trade-off problem in Sobol' index estimation for stochastic models and they allow to show that the strategy proposed in this paper performs well.

Conclusion

To balance the tradeoff between repetitions and explorations in the estimation of Sobol' indices for stochastic models, a regularized estimator was devised and its mean-squared error bounded, approximated and minimized. A connection to functions of nested Monte Carlo estimators was made, for which generic results were obtained, especially useful when the Hessian matrix is unbounded. To estimate Sobol' indices in practice, an algorithm that adpats to the intrinsic randomness of the stochastic model was designed and illustrated on simulations, where the importance of a good balance between explorations and repetitions was observed. In the theoretical results, the regularization parameter, though arbitrary, has to be fixed. In the future, it would be of primary interest to let it vanish as the number of explorations and repetitions go to infinity, allowing one to get a genuine convergence rate for the mean-squared error of the regularized estimator. A further step would be to get a minimax low bound. Also, it would be interesting to compare the results with multilevel Monte Carlo methods [START_REF] Mycek | Multilevel Monte Carlo covariance estimation for the computation of Sobol' Indices[END_REF]; [START_REF] Giles | Multilevel nested simulation for efficient risk estimation[END_REF]). Furthermore, it could be interesting to couple the iterative estimation approach of [START_REF] Gilquin | Iterative estimation of Sobol' indices based on replicated designs[END_REF] to the algorithms implemented in this study in order to build an adaptive version which could perform estimation with respect to a given precision. Another estimation approach worth considering is that based on Chatterjee's empirical correlation coefficient [START_REF] Chatterjee | A New Coefficient of Correlation[END_REF][START_REF] Gamboa | Global sensitivity analysis: A novel generation of mighty estimators based on rank statistics[END_REF].

Appendix A Proof of Proposition 1

Using convexity inequality, for all h ∈ (0, 1) and m ≥ 1, it holds:

E g h θ -g (θ) 2 ≤ 2E g h ( θ) -g h (µ m ) 2 + 2 (g h (µ m ) -g (θ)) 2 .
Applying a Taylor-Lagrange expansion to g h at points θ and µ m yields:

E g h ( θ) -g h (µ m ) 2 = E ∇g h (µ m ) ⊤ θ -µ m + 1 2 θ -µ m ⊤ ∇ 2 g h λ θ + (1 -λ)µ m θ -µ m 2 ≤ 2E ∇g h (µ m ) ⊤ θ -µ m 2 + 1 2 E θ -µ m ⊤ ∇ 2 g h λ θ + (1 -λ)µ m θ -µ m 2 ,
for some λ ∈ (0, 1). Thus:

E g h θ -g h (µ m ) 2 E ∇g h (µ m ) ⊤ θ -µ m 2 ≤ 2 + 1 2 E θ -µ m ⊤ ∇ 2 g h λ θ + (1 -λ)µ m θ -µ m 2 E ∇g h (µ m ) ⊤ θ -µ m 2 Let p n,m (h) := sup λ∈[0,1] E θ -µ m ⊤ ∇ 2 g h λ θ + (1 -λ)µ m θ -µ m 2 E ∇g h (µ m ) ⊤ θ -µ m 2 , (16)
then

E g h ( θ) -g h (µ m ) 2 ≤ (2 + 1 2 p n,m (h))V n,m (h) ≤ 2(1 + p n,m (h))V n,m (h).
Now, let us show that p n,m (h) → 0 as n, m → ∞. Using Cauchy-Schwarz inequality yields that:

p n,m (h) ≤ sup λ∈[0,1] E ∇ 2 g h λ θ + (1 -λ)µ m 4 F × E θ -µ m 8 E ∇g h (µ m ) ⊤ θ -µ m 2 .
As n, m → ∞, the first term of the product above is bounded by a constant independent of h, uniformly in λ, by Assumption 2. The second term is given by:

E θ -µ m 8 E ∇g (µ m + hu) ⊤ θ -µ m 2 (17) = E θ -µ m 8 E ∇g (µ m ) ⊤ θ -µ m 2 × E ∇g (µ m ) ⊤ θ -µ m 2 E ∇g (µ m + hu) ⊤ θ -µ m 2 ≤ E θ -µ m 8 E ∇g (µ m ) ⊤ θ -µ m 2 × sup h∈(0,1) E ∇g (µ m ) ⊤ θ -µ m 2 E ∇g (µ m + hu) ⊤ θ -µ m 2 . ( 18
)
Lemma A.1.

E θ -µ m 8 E ∇g (µ m ) ⊤ θ -µ m 2 = o(1), (19)
as n, m → +∞

The first factor in the right-hand side in (18) does not depend on h and is of order o(1) as stated in Lemma A.1. Moreover, since ∇g is continuous in parameter h, we have that the second factor in the right-hand side of (18) is O(1) as n, m → ∞. Therefore, p n,m (h) = o(1), and hence

E g h ( θ) -g(θ) 2 ≤ 4 (1 + p n,m (h)) V n,m (h) + B m (h) 2 ,
where B m (h) = g h (µ)g(θ) and lim n,m→+∞ p n,m (h) = 0.

Proof of Lemma A.1

First, let us bound the numerator of the ratio in Equation ( 19); we have

E θ -µ m 8 ≤ 27 3 j=1 E θ j -µ mj 8 ,
where θ j and µ mj denote the jth component of θ and µ m , respectively. By [START_REF] Marcinkiewicz | Sur les fonctions indépendantes[END_REF] and Jensen inequalities, we have for every j that

E θ j -µ mj 8 ≤ B 8 n 4 E θ (1) mj -µ mj 8 ,
where here B 8 is a universal constant. The case j = 2 is the simplest. Notice that µ m2 does not depend on m, then the expansion of E(| Q m (X (1) )-µ m2 | 8 ) through Newton formula yields terms of the form µ k m2 E( Q 8-k m ), k = 0, . . . , 8. Using Lemma C.1 in Section C provides that those terms are polynomial in m -1 .

Let us deal with the case j = 1. Expanding the power 8 through Newton's formula and bounding its terms yields

E Q m (X (1) ) 2 -µ m1 8 ≤ (µ 8 m1 ∨ 1) 8 4 E Q m (X (1) ) 16 + 1 . ( 20
)
Denoting ϕ(X (1) , Z (1,k) ) = Y (1,k) , we have k16) .

E Q m (X (1) ) 16 = E   1 m m k=1 Y (1,k) 16   = 1 m 16 m k1,...,k16=1 E Y (1,k1) • • • Y (1,
The expectation in the right-hand side is symmetric in k 1 , . . . , k 16 , and hence, from Lemma C.1, the sum is a polynomial in m of degree 16. Therefore, the right-hand side in ( 20) is bounded uniformly in m.

Let us deal with the case j = 3. Proceeding as in (20), we have

E Q m (X (1) ) Q m ( X (1) ∼u , X (1) u ) -µ m3 8 ≤ (µ 8 m3 ∨ 1) 8 4 E Q m (X (1) ) Q m ( X (1) ∼u , X (1) u ) 8 + 1 ≤ (µ 8 m3 ∨ 1) 8 4 E 1 2 Q m (X (1) ) 16 + 1 2 Q m ( X (1) ∼u , X (1) u ) 16 + 1 ,
and this is also bounded uniformly in m. (Again by Lemma C.1.)

We now deal with the root of the denominator of the ratio in Equation ( 19). We have

E ∇g (µ m ) T θ -µ m 2 = 3 j1,j2=1 ∇g(µ m ) j1 ∇g(µ m ) j2 E( θ -µ m ) j1 ( θ -µ m ) j2 = 1 n 3 j1,j2=1 ∇g(µ m ) j1 ∇g(µ m ) j2 E θ (1) mj1 θ (1) mj2 -µ mj1 µ mj2 . (21)
The infimum of the sum in ( 21) is reached for some m and greater than zero. Therefore, the numerator in Equation ( 19) is less than 1/n 4 times a constant not depending on m or n and the denominator is equal to 1/n 2 times a quantity greater than zero. Therefore, the supremum over m of the ratio in Equation ( 19) is of order O(n -2 ). The proof is complete.

B Proof of Theorem 1

The following lemma will be needed:

Lemma B.1. For all α > 0 and all h ∈ (0, 1),

lim n,m→+∞ E   1 h + θ 1 -θ 2 2 α   = 1 (h + Var (Q(X))) α ≤ Var (Q(X)) -α . (22)
Note that the function g is infinitely differentiable over its domain D and then its gradient is given by:

∇g (θ 1 , θ 2 , θ 3 ) = - θ 3 -θ 2 2 (θ 1 -θ 2 2 ) 2 2θ 2 (θ 3 -θ 1 ) (θ 1 -θ 2 2 ) 2 , 1 θ 1 -θ 2 2 ⊤ .
Furthermore, the hessian matrix of g yields:

∇ 2 g(θ 1 , θ 2 , θ 3 ) =     2(θ3-θ 2 2 ) (θ1-θ 2 2 ) 3 2θ2(θ1-2θ3+θ 2 2 ) (θ1-θ 2 2 ) 3 -1 (θ1-θ 2 2 ) 2 2θ2(θ1-2θ3+θ 2 2 ) (θ1-θ 2 2 ) 3 2(θ3-θ1)(θ1+3θ 2 2 ) (θ1-θ 2 2 ) 3 2θ2 (θ1-θ 2 2 ) 2 -1 (θ1-θ 2 2 ) 2 2θ2 (θ1-θ 2 2 ) 2 0     . For any (θ 1 , θ 2 , θ 3 ) ∈ D, the matrix ∇ 2 g(θ 1 , θ 2 , θ 3 ) is under the form ∇ 2 g(θ 1 , θ 2 , θ 3 ) = B(θ 1 , θ 2 , θ 3 )/(θ 1 -θ 2 2 ) 3 where B(θ 1 , θ 2 , θ 3 ) is the matrix: B(θ 1 , θ 2 , θ 3 ) =   2(θ 3 -θ 2 2 ) 2θ 2 (θ 1 -2θ 3 + θ 2 2 ) -(θ 1 -θ 2 2 ) 2θ 2 (θ 1 -2θ 3 + θ 2 2 ) 2(θ 3 -θ 1 )(θ 1 + 3θ 2 2 ) 2θ 2 (θ 1 -θ 2 2 ) -(θ 1 -θ 2 2 ) 2θ 2 (θ 1 -θ 2 2 ) 0   .
Notice that B(θ 1 , θ 2 , θ 3 ) includes only multivariate polynomials of variables θ 1 , θ 2 and θ 3 . Let us check Assumption 2. Let λ ∈ [0, 1] and h ∈ (0, 1). We have

∇ 2 g λ θ + (1 -λ)µ m + hu = B λ θ + (1 -λ)µ m + hu h + λ θ 1 + (1 -λ)µ m1 -λ θ 2 + (1 -λ)µ m2 2 3 .
Thus:

E sup λ ∇ 2 g λ θ + (1 -λ)µ m + hu 4 F = E      sup λ B λ θ + (1 -λ)µ m + hu 4 F h + λ θ 1 + (1 -λ)µ m1 -λ θ 2 + (1 -λ)µ m2 2 12      ≤ E      sup λ 1 h + λ θ 1 + (1 -λ)µ m1 -λ θ 2 + (1 -λ)µ m2 2 24      × E sup λ B λ θ + (1 -λ)µ m + hu 8 F ≤ E   sup λ 1 h + λ θ 1 -θ 2 2 + (1 -λ)(µ m1 -µ 2 m1 ) 24    (by convexity inequality) × E sup λ B λ θ + (1 -λ)µ m + hu 8 F ≤ E   sup λ λ h + θ 1 -θ 2 2 24    + sup λ 1 -λ (h + µ m1 -µ 2 m1 ) 24
(by convexity inequality)

× E sup λ B λ θ + (1 -λ)µ m + hu 8 F ≤ E    1 h + θ 1 -θ 2 2 24    + 1 (h + µ m1 -µ 2 m1 ) 24 × E sup λ B λ θ + (1 -λ)µ m + hu 8 F ≤ E    1 h + θ 1 -θ 2 2 24    + sup h∈(0,1) 1 (h + µ m1 -µ 2 m1 ) 24 × sup h∈(0,1) E sup λ B λ θ + (1 -λ)µ m + hu 8 F .
One should remark that sup h∈(0,1)

1 (h+µm1-µ 2 m1 ) 24 ≤ 1 (Var(E[ϕ(X,Z)|X])) 24 < +∞.
Moreover, the matrix B is composed with polynomials of three variables. Since E Q(X) 16 < +∞ then by using Lemma C.1 and by continuity of polynomial functions, it yields that sup h∈(0,1) E sup λ∈[0,1] B λ θ + (1λ)µ m + hu 8 F is bounded. Finally, by relying on Lemma B.1, E(h + θ 1θ 2 2 ) -24 is a bounded by 1 (Var(E[ϕ(X,Z)|X])) 24 as n, m → +∞. Therefore, Assumption 2 is satisfied.

Proof of Lemma B.1

Let h ∈ (0, 1) be fixed. The function β h : x → 1/(h + x) α is continuously differentiable such that its first derivative is uniformly bounded on R + by 1/h α+1 then it is Lipschitz. Therefore: Straightforwardly:

E β h ( θ 1 -θ 2 2 ) -β h (θ 1 -θ 2 2 ) 2 ≤ 1 h 2α+2 E θ 1 -θ 2 2 -θ 1 + θ 2
lim n,m→∞ E 1 (h + ( θ 1 -θ 2 2 ) α = 1 (h + θ 1 -θ 2 2 ) α ≤ 1 Var(Q(X)) α .

C A lemma

Lemma C.1. Let X (1) , • • • , X (n) be n i.i.d. copies of X and Z (1,1) , • • • , Z (n,m) be n × m i.i.d. copies of Z such that (X (1) , • • • , X (n) ) and (Z (1,1) , • • • , Z (n,m) ) are independent. Then, for all q ∈ N: m -q E m k=1 ϕ(X (1) , Z (1,k) )

q is polynomial in m -1 of degree q -1 with constant E (E [ϕ(X, Z) | X]) q .

It holds that: 1) , Z (1,kq) ) 1) , Z (1,kq) ) .

E 1 m m k=1 ϕ(X (1) , Z (1,k) ) q = 1 m q E   m k1=1 • • • m kq=1 ϕ(X (1) , Z (1,k1) ) • • • ϕ(X (
  = 1 m q m k1=1 • • • m kq=1 E ϕ(X (1) , Z (1,k1) ) • • • ϕ(X (
Denote by λ : {1, . . . , m} q → N the map which with each k := (k 1 , . . . , k q ) associates the number of distinct indices among k 1 , . . . , k q . If 1 ≤ l ≤ q then denote by ρ l : λ -1 (l) → {1, . . . , q} l the map which with each k ∈ λ -1 (l) associates (r 1 , . . . , r l ), where r i = |{j : k j = k ji }| for every i = 1, . . . , l and k j1 , . . . , k j l are the distinct indices found among k 1 , . . . , k q . Obviously, r 1 + • • • + r l = q. We have Now, since E ϕ(X (1) , Z (1,k1) ) • • • ϕ(X (1) , Z (1,kq) ) = E ϕ(X (1) , Z (1,kj 1 ) ) r1 • • • ϕ(X (1) , Z (1,kj l ) ) r l = E l s=1 E ϕ(X (1) , Z (1,kj s ) ) rs | X (1) is symmetric in r 1 , . . . , r l , it holds that E ϕ(X (1) , Z (1,kj s ) ) rs | X (1) where c(l, (r 1 , . . . , r l ), m) = q r 1 qr 1 r 2

• • • q -r 1 -• • • -r l-1 r l m(m -1) • • • (m -l + 1). ( 24
)
Notice that the expression in the right-hand side of ( 24) is invariant by permutation of r 1 , . . . , r l . Therefore, the sum ( 23) is a polynomial in m of degree q with constant zero and hence E 1 m m k=1 ϕ(X (1) , Z (1,k ) q is a polynomial in 1 m of degree q -1 with constant lim m→+∞ E 1 m m k=1 ϕ(X (1) , Z (1,k ) q = E (E [ϕ(X, Z) | X]) q .

D The Marcinkiewicz-Zygmund inequality Theorem D.1 [START_REF] Marcinkiewicz | Sur les fonctions indépendantes[END_REF]). Let U 1 , • • • , U n be i.i.d. random variables such that E(U 1 ) = 0 and E|U 1 | q < +∞, where 1 ≤ q < +∞. There exist A q and B q depending only on q such that:

A q E   n i=1 |U i | 2 q 2   ≤ E n i=1 U i q ≤ B q E   n i=1 |U i | 2 q 2  
Furthermore, there exists C q independent from n such that:

E 1 n n i=1 U i q ≤ C q n q 2 . ( 25
)

Figure 1 :Figure 2 :

 12 Figure1: Boxplots of sensitivity index estimates of X 1 and X 2 with respect to the linear model and for different values of T . The standard deviation is σ = 1 and regularization parameter h = 10 -2 . Three strategies of choice of m are compared: m = 5 (in red), m = m opt given by the trade-off strategy of Algorithm 2 (in green) and m = T 1/2 (in blue). The red dashed lines showed the true sensitivity index values, in particular, in this setting S 1 = 0.2 and S 2 = 0.8

Figure 3 :Figure 4 :

 34 Figure 3: Boxplots of sensitivity index estimates of X 1 , X 2 and X 3 with respect to the stochastic version of Ishigami function (with b = 0.05) for different values of T . Three strategies of choice of m are compared: m = 5 (in red), m = m opt given by the trade-off strategy of Algorithm 2 (in green) and m = T 1/2 (in blue).

2 ,

 2 Var( θ 1 ) + Var( θ 2 2 ) + E( θ 2 2 )θ 2 2 using convexity inequality. Based on Marcinkiewicz-Zygmund inequality (see Theorem D.1) and Lemma C.1, it follows that lim n,m→∞ (Var( θ 1 ) + Var( θ 2 2 ) + E( θ 2 2 )

  X (1) , Z (1,k1) ) • • • ϕ(X (1) , Z (1,kq) ...,r l )∈{1,...,q} l :r1+•••+r l =q   k∈λ -1 (l):ρ l (k)=(r1,...,r l )

(

  r1,...,r l )∈{1,...,q} l :r1+•••+r l =q   k∈λ -1 (l):ρ l (k)=(r1,...,r l ) f (k)   = c(l, (r 1 , . . . , r l ), m)E l s=1

Acknowledgment

We thank a anonymous reviewer for their helpful comments. We thank Clémentine Prieur for reading the first version of this paper and for the pertinent remarks and references she suggested for its improvement.

Algorithm 1 Completing model evaluations


T P Tuyen Dao

Lubomir Vezenkov

Gil Subra

Muriel Amblard

Vincent Ladmiral

Mona Semsarilar

Morphological evolution of poly(glycerol monomethacrylate)-stat-glycine-phenylalaninephenylalanine-methacrylamide-b-poly(2hydroxypropylmethacrylate) (P(GMA-stat-(MAm-GFF))-b-PHPMA) block copolymer nanoobjects via polymerization-induced self-assembly †

Here we report a study on the reversible addition-fragmentation chain transfer (RAFT) polymerizationinduced self-assembly (PISA) of self-assembling peptide-containing diblock copolymers. The self-assembling tripeptide (SAP) monomer used (MAm-GFF) is a methacrylamide-functionalized Gly-Phe-Phe sequence. The incorporation of this self-assembling GFF tripeptide into the stabilizing PGMA segment of a well-known poly(glycerol monomethacrylate)-b-poly(2-hydroxypropyl methacrylate) block copolymer (PGMA-b-PHPMA) has been shown in a previous study to induce the formation of complex selfassembled structures. The present article explores the effect of varying the block length (both hydrophilic and hydrophobic blocks), block composition (different SAP units in the hydrophilic block), solvent composition and temperature on the morphologies of PISA synthesized P(GMA-stat-(MAm-GFF))-b-PHPMA.

The polymerization of HPMA from the SAP-containing P(GMA-stat-(MAm-GFF)) stabilizing block led to well-defined diblock copolymers, which self-assembled into diverse morphologies ranging from vesicles, worms, or micelles to unreported highly anisotropic structures. Moreover, these morphologies were found to be thermoresponsive. At 70 °C, more classic morphologies were observed as hydrogen bonds are weakened at high temperatures. However, these structures transform into unprecedented morphologies (urchin-like) once cooled to 4 °C probably due to the reduced hydrophilicity of the PHPMA block.

Introduction

Self-assembly of amphiphilic molecules in solution is a very promising bottom-up approach in nanomaterials science. It has found many potential applications in the fields of photonics, electronics, sustainable energy, bioengineering or drug delivery for instance. [1][2][3][START_REF] Zhang | Self-assembly of photonic crystals from polymer colloids[END_REF][START_REF] Liu | Preparation and Characterization of Nanomaterials for Sustainable Energy Production[END_REF][START_REF] Yadav | Nanoscale Self-Assembly for Therapeutic Delivery[END_REF][START_REF] Stupp | Self-assembly of biomolecular soft matter[END_REF][START_REF] Shimizu | Self-organized nanotube materials and their application in bioengineering[END_REF] Besides the well-known morphologies (micelles, worm-like micelles or vesicles), the discovery of more complex architectures (such as tubes, rods, platelets, sheets, wires, ribbons, interwound fibrils, bundles, etc.) has stimulated many studies aiming to understand and control their formation. Molecular self-assembly results from the combination of various weak forces such as hydrophobic interactions, electrostatic interactions, hydrogen bonding and π-π stacking.

Block copolymer self-assembly has usually been performed using the post-polymerization approach whereby preformed block copolymers are self-assembled in a selected solvent. This approach is, however, limited to relatively low concentrations (typically less than 1%), and is both time-and labor-intensive as it involves multiple steps ( polymerization, purification, and self-assembly), is not easily reproducible and often does not lead to a unique type of morphology. Polymerization-induced self-assembly (PISA), proceeding via the simultaneous formation and self-assembly of amphiphilic block copolymers, is a powerful method to prepare various morphologies at high concentrations (up to 50% solid content), fast (the polymeriz- † Electronic supplementary information (ESI) available. See DOI: https://doi.org/ 10.1039/d2py01203k a Institut Européen des Membranes, IEM, Univ. Montpellier, CNRS, ENSCM, Montpellier, France. E-mail: mona.semsarilar@umontpellier.fr b Institut Charles Gerhardt Montpellier, ICGM, Univ. Montpellier, CNRS, ENSCM, Montpellier, France ation time) and with high reproducibility. [START_REF] Penfold | Emerging Trends in Polymerization-Induced Self-Assembly[END_REF][START_REF] Liu | Polymerization techniques in polymerization-induced self-assembly (PISA)[END_REF][11][START_REF] Canning | A Critical Appraisal of RAFT-Mediated Polymerization-Induced Self-Assembly[END_REF][START_REF] D'agosto | RAFT-Mediated Polymerization-Induced Self-Assembly[END_REF] In addition, PISA can be used with different types of polymerization techniques and in various media. [START_REF] Derry | Polymerizationinduced self-assembly of block copolymer nanoparticles via RAFT non-aqueous dispersion polymerization[END_REF][START_REF] Warren | Polymerization-Induced Self-Assembly of Block Copolymer Nano-objects via RAFT Aqueous Dispersion Polymerization[END_REF][START_REF] Derry | Industrially-relevant polymerization-induced self-assembly formulations in non-polar solvents: RAFT dispersion polymerization of benzyl methacrylate[END_REF] Few studies have reported the synthesis of self-assembled nanostructures based on peptides via PISA. [START_REF] Ladmiral | Synthesis and characterization of poly(amino acid methacrylate)-stabilized diblock copolymer nano-objects[END_REF][START_REF] Luppi | Antimicrobial polylysine decorated nano-structures prepared through polymerization induced self-assembly (PISA)[END_REF][START_REF] Sun | Proapoptotic Peptide Brush Polymer Nanoparticles via Photoinitiated Polymerization-Induced Self-Assembly[END_REF] Peptides are able to self-assemble via a range of non-covalent interactions. This ability is responsible for the secondary structure of proteins. We recently studied the effect of the insertion of selfassembling peptide (SAP) moieties into diblock copolymers on their PISA behaviour. Our group started this approach by inserting SAPs into very well-studied PGMA-b-PHPMA amphiphilic copolymers to investigate their influence on the final PISA-generated structures. 1 This study showed that the insertion of only a few monomer units containing the SAP (tripeptide of Gly-Phe-Phe abbreviated as GFF) could result in the formation of long anisotropic fibrous morphologies, confirming the decisive influence of the GFF SAP on the self-assembly of P(GMA-stat-(MAm-GFF))-b-PHPMA.

In another report, we used SAP containing monomers as the core-forming blocks, either on their own (forming a pure polypeptide block) or as a statistical copolymer with HPMA and GMA. [START_REF] Dao | Nano-assemblies with core-forming hydrophobic polypeptide via polymerization-induced self-assembly (PISA)[END_REF] In all cases, placing the SAP-containing monomer in the core resulted in the formation of a range of interesting morphologies (fibre-, flake-, and leaf-like) never reported previously.

Here, we report a more detailed study of the P(GMA-stat-GFF)-b-PHPMA formulation. Our previous report dealt with only one stabilizing block composition. 1 Here, different parameters influencing the morphology and its evolution, namely, the number of SAP monomers in the CTA, the block length, solvent ( pure water and a mixture of water/alcohol) and temperature were studied. The synthetic approach in this work is illustrated in Scheme 1.

Experimental section

Materials

Glycerol monomethacrylate (GMA), 2-hydroxypropyl methacrylate (HPMA), 4,4′-azobis(4-cyanopentanoic acid) (ACVA) and all organic solvents were purchased from Sigma-Aldrich. DMSO-d 6 was purchased from Eurisotop. MAm-GFF and 4-cyano-4-(2phenylethanesulfanylthiocarbonyl)-sulfanyl pentanoic acid (PETTC) were prepared according to our previous reports. 1,[START_REF] Semsarilar | Anionic Polyelectrolyte-Stabilized Nanoparticles via RAFT Aqueous Dispersion Polymerization[END_REF] Syntheses of P(GMA-stat-MAm-GFF) macromolecular chain transfer agents P(GMA x -stat-(MAm-GFF) y ) macromolecular chain transfer agents (mCTAs) were synthesized via RAFT solution copolymerization according to the following procedure. A round bottom Schlenk flask equipped with a stirrer bar was charged with GMA (x eq.), MAm-GFF (y eq.), PETTC RAFT agent (1 eq.), ACVA initiator (0.1 eq.) and DMF at the volume needed to solubilize all of the solids (the final concentration of solids was about 20% w/w). The flask was subsequently sealed with a rubber septum, placed in an ice bath and degassed by bubbling nitrogen for 30 min. The flask was then immersed in a preheated oil bath at 70 °C under stirring for 24 h. The copolymerization was quenched by opening the flask to air. GMA and MAm-GFF conversions were determined from the 1 H NMR spectrum in DMSO-d 6 (Fig. S1 †). The reaction solution was diluted in DMF and precipitated in cold THF. The yellow precipitate was then redispersed in water and freeze-dried. The final product was a light yellow powder. The mean DP was determined by 1 H NMR and the molar mass distribution was determined using SEC-HPLC in DMAc or DMF.

Preparation of diblock P(GMA-stat-(MAm-GFF))-b-PHPMA via dispersion RAFT-PISA

A representative PISA synthesis via RAFT dispersion polymerization of HPMA with the P(GMA-stat-(MAm-GFF)) macro-CTA was carried out as described below. For a target composition of P(GMA 24 -stat-(MAm-GFF) 3 )-b-PHPMA 60 , the purified P(GMA 24 -stat-(MAm-GFF) 3 ) macro-CTA (100 mg, 0.018 mmol), HPMA (155.5 mg, 1.08 mmol) and ACVA initiator (1.0 mg, 0.0036 mmol) were weighed into a 10 mL glass vial containing a stirrer bar. 2.3 mL Milli-Q water was then added to obtain a reaction solution at 10% w/w. The vial was sealed, degassed for 30 minutes and placed in an oil bath previously set at 70 °C using a stirring hotplate. After 24 h, the reaction was opened to air, and cooled down to ambient temperature. The HPMA conversion was always 100% with no remaining HPMAvinyl signals in 1 H NMR spectra. Other aliquots were taken to perform SEC, DLS and TEM.

In the following text, G x M y denotes the P(GMA x -stat-(MAm-GFF) y ) macro CTAs and G x M y H z denotes the P(GMA x -stat-(MAm-GFF) y )-b-PHPMA z formulation.

Characterization

Proton nuclear magnetic resonance spectroscopy ( 1 H NMR). 1 H NMR spectra were recorded using a 300 MHz Bruker Avance-300 spectrometer, processed and analyzed with MestReNova 9.0 software.

Size exclusion chromatography (SEC-HPLC). Polymers were analysed using a Varian 390-LC system equipped with 2 PL1113-6300 ResiPore 300 × 7.5 columns and a 390LC PL0390-0601 refractive index detector (RI). The mobile phase was DMF with 0.1% w/w LiBr adjusted at a flow rate of 1 mL min -1 while the columns were thermostated at 70 °C. Calibration was performed with near-monodisperse poly(methyl methacrylate) (PMMA) standards ranging from 550 to 1 568 000 g mol -1 (EasiVial-Agilent).

Dynamic light scattering (DLS). The hydrodynamic radii R H were analysed by dynamic light scattering at 90°using a Litesizer TM 500 Anton Paar. Samples were prepared by diluting directly from PISA-suspension 100-fold or dispersing the resulting polymers in filtered Milli-Q water at 0.1% w/w. The measurements were performed at 20 °C or at a series of temperatures from 20 °C to 70 °C with 10 °C intervals, and 30 min of equilibration was set at the target temperature prior to each measurement.

Transmission electron microscopy (TEM). TEM images were acquired using either a JEOL 1200 EXII-120 kV or JEOL 1400 P+ -120 kV microscope. To prepare the TEM grid, 10 µL of 100-fold diluted sample from the PISA suspension (0.1% w/w) was deposited onto the grid for 60 s and then blotted with filter paper to remove excess solution. Afterward, the sampleloaded grid was stained with 7 µL of 1% ammonium molybdate solution for 20 s. Excess stain was also removed with filter paper. The grid was allowed to dry under a hood for 5 minutes.

TEM grids for imaging at 70 °C or 4 °C were prepared as follows: a 10% w/w PISA suspension was incubated at the desired temperature for 1 h, diluted to 0.1% w/w at the same temperature and this solution was deposited on a TEM grid and stained as described above. Adequately thermostated glassware was used for sample preparation.

Results and discussion

P(GMA-stat-MAm-GFF) macromolecular chain transfer agents (mCTAs)

Copolymers of the SAP-containing monomers MAm-GFF and GMA of various compositions were synthesized by RAFT in DMF using PETTC as the RAFT agent and ACVA as the initiator. The solid content in the reaction was about 20 wt% due to the low solubility of MAm-GFF in DMF. The copolymerization was carried out for 24 h. This duration led to the full conversion of GMA while 80% of MAm-GFF was copolymerized. According to the kinetics study performed using in situ 1 H NMR and reported in our previous work, 1 MAm-GFF was indeed consumed slower than GMA. However, a statistical P(GMA-stat-(MAm-GFF)) copolymer was formed.

The resulting mCTAs (mCTA 1, mCTA 2 and mCTA 3) were characterized by 1 H NMR and SEC (Fig. S1, S2 † and Table 1). The copolymerization proceeded as expected with a satisfactory level of control and the resulting mCTA had relatively low dispersity (<1.4). These 3 mCTAs afforded a clear solution in Milli-Q water at 10 wt%, enabling the aqueous dispersion PISA protocols.

TEM and DLS analyses were also performed to investigate the mCTA aqueous solutions used for the PISA experiments. 0.1 wt% DLS and TEM samples were prepared from 10 wt% stock solutions in Milli-Q water. As shown in our previous work, mCTA 2 formed large bundles of fibers/ribbons with diameters of 60-500 nm and an average length of 10 µm (Fig. 1). mCTA 1 also formed large bundles of fibers/ribbons. However, as shown in Fig. 1, the structures are not exactly identical. These structures likely still exist at 70 °C based on our previous study 1 and the GMA/MAm-GFF ratios of mCTA 1 and mCTA 2 are similar. In contrast, no clear structure was found for the mCTA 3 aqueous solution. This suggests that this mCTA was molecularly dissolved at 0.1 wt% as expected given the GMA/ MAm-GFF ratio (200/9) of this copolymer. Light scattering data confirmed the presence of objects with very high dispersity in both mCTA 1 and 2 solutions (Fig. S3 †) in agreement with the anisotropic morphologies observed by TEM. No morphology was observed for the longer macro-CTA with fewer MAm-GFF units (mCTA 3), indicating full solubility under the tested conditions.

RAFT aqueous dispersion polymerisation of HPMA using the MAm-GFF-containing mCTAs Amphiphilic P(GMA x -stat-(MAm-GFF) y )-b-PHPMA z (G x M y H z ) block copolymers (BCPs) featuring different PHPMA block lengths were synthesized from the three MAm-GFF-containing mCTAs under aqueous RAFT PISA conditions. The polymerization proceeded without any issues, as expected, reaching full conversion (>99%) within 24 h and producing relatively welldefined block copolymers (Đ < 1.35, Fig. 2). The M n and Đ values of all diblock copolymers are given in Table S1. † a G x M y : x and y correspond to the DP of GMA and MAm-GFF calculated from the 1 H NMR spectra (Fig. S1 †), respectively. In all cases, the polymerization resulted in stable suspensions when low to moderate DP PHPMA s were targeted. These BCP suspensions were stable for months (>12). However, the block copolymers with a higher DP PHPMA phase separated over time into an excluded liquid phase and a viscous phase that could not be re-dispersed even after continuous vigorous stirring for up to 7 days (Fig. 3). Among the BCPs prepared from mCTA 1, only G 24 M 3 H 14 formed a stable suspension. G 24 M 3 H 29 , G 24 M 3 H 58 and G 24 M 3 H 108 led to phase separation within two hours after polymerization. Similarly, the suspensions of G 65 M 7 H 28 and G 65 M 7 H 42 were stable while those of G 65 M 7 H 54 and G 65 M 7 H 63 phase-separated one hour after the polymerization, and G 65 M 7 H 200 formed a gel during the PISA process (at 70 °C). In contrast, the diblock copolymers prepared using mCTA 3 formed stable suspensions. The cause of the phase separation is certainly due to the presence of SAP moieties, since such a phase separation has never been reported for any of the PGMA-b-PHPMA diblock copolymers synthesized under the same conditions (aqueous PISA at 10 wt%) even when very long PHPMA blocks (PGMA 47 -b-PHPMA 310 , PGMA 78 -b-PHPMA 500 , and PGMA 112 -b-PHPMA 800 ) [START_REF] Warren | Polymerization-Induced Self-Assembly of Block Copolymer Nano-objects via RAFT Aqueous Dispersion Polymerization[END_REF][START_REF] Warren | Critical Dependence of Molecular Weight on Thermoresponsive Behavior of Diblock Copolymer Worm Gels in Aqueous Solution[END_REF] were targeted. As reported previously the presence of only a few units (3-7) of MAm-GFF is enough to profoundly change the morphologies of the self-assemblies of PGMA-b-PHPMA BCPs and the stabilities of their suspensions. Without doubt, the colloidal instability of these nanoparticle solutions is due to the reduced solubility of the SAP containing mCTAs. When the length of the mCTA is increased and the ratio of MAm-GFF/GMA is reduced as in the case of mCTA 3, stable colloidal solutions could be prepared even at very high DPs of PHPMA (≥500).

TEM analyses were performed for all diblock copolymers by simply diluting the stable suspension 100-fold or by re-dispersing the phase-separated solids (viscous phase that sedimented at the bottom of the vials) in Milli-Q water at 1 mg mL -1 . Our initial study of this system had shown that G 65 M 7 H 28 BCP PISA resulted in anisotropic rigid fibres that formed a fibrous network upon drying. 1 Here, as shown in Fig. 4, the selfassembled morphologies observed post-mortem (at 25 °C) for the different G x M y H z series are diverse. In the G 24 M 3 H z series, only G 24 M 3 H 14 formed colloidally stable solutions. This series generally self-assembled into fibers of about 5 µm length, independently of their final colloidal stability. Increasing the length of the PHPMA block to 29 resulted in the formation of clustered fibers that formed separated bundles when z was increased to 58 and 108 with significantly more branching. This change in the morphology is no surprise as increasing the hydrophobicity of the BCP often results in the formation of higher order morphologies. However, maintaining the bundled fibrous networks (established as the specific feature of PISA formulations containing the MAm-GFF units) up to very high DPs of hydrophobic PHPMA (only 3 units of MAm-GFF with 108 units of HPMA) is remarkable. As reported in our previous work, 1 in the absence of MAm-GFF, an increase in the length of PHPMA resulted in morphological changes towards worms and vesicles.

For the G 65 M 7 H z series (z ranging from 42 to 63), once again, bundles of fibers were observed (Fig. 4). The individual fibers seem identical in all the samples but they tended to arrange into different spherulitic structures (network-like for z = 42, bundles with a darker centre for z = 54 and sea urchinlike for z = 63). Finally, 55 nm diameter (measured by DLS) spheres were obtained for G 65 M 7 H 200 , suggesting almost complete suppression of the effect of the SAP since in the absence of the peptide sequence (G 65 H 200 ), spheres of 87 nm were observed. [START_REF] Li | RAFT Synthesis of Sterically Stabilized Methacrylic Nanolatexes and Vesicles by Aqueous Dispersion Polymerization[END_REF] It is important to note that some of these G 65 M 7 H x formulations did not form stable suspensions. The TEM images of these samples were obtained by re-dispersing the colloidally unstable solution in water at a lower concentration (0.1 wt%) when possible, or by sampling and diluting the viscous/solid phase that sedimented at the bottom of the vials in Milli-Q water (at 1 mg mL -1 ).

The G 200 M 9 H z series prepared from water-soluble mCTA 3 behaved quite differently. Only G 200 M 9 H 51 led to spherulitic structures (Fig. 4). The other BCPs, G 200 M 9 H 102 , G 200 M 9 H 250 , and G 200 M 9 H 510 , self-assembled into more classical morphologies: fused spheres (dimers and trimers), spheres and vesicles, respectively (Fig. 4). These are the typical morphologies observed for the PGMA-b-PHPMA BCP suggesting that the self-assembly of G 200 M 9 H x with x = 102-510 was not directed by the interactions of the GFF moieties but by the hydrophobicity/hydrophilicity balance as in the case of the PGMA-b-PHPMA system. These results underline once again that the presence of very few units of MAm-GFF can have a significant impact on the morphology of PGMA-b-PHPMA polymerization-induced self-assembled objects. This impact of the SAP is all the more pronounced when the PGMA block is shorter. However, when the length of the stabilising block is relatively long and the ratio of GMA/MAm-GFF is low the effect of the peptide sequence is significantly reduced or completely suppressed particularly at high DPs of the core-forming block (PHPMA).

Temperature-induced morphological changes

Besides the composition of the block copolymers, experimental conditions such as ionic strength, pH and temperature are key parameters to access different self-assembled morphologies. [START_REF] Zhang | pH-responsive self-assembly of fluorophore-ended homopolymers[END_REF][START_REF] Wakabayashi | pH-Responsive Self-Assembly of Designer Aromatic Peptide Amphiphiles and Enzymatic Post-Modification of Assembled Structures[END_REF][START_REF] Bennevault | Temperature sensitive supramolecular self assembly of per-6-PEO-β-cyclodextrin and α,ω-di-(adamantylethyl)poly(N-isopropylacrylamide) in water[END_REF][START_REF] Zhou | Polymerization induced self-assembly: tuning of morphology using ionic strength and pH[END_REF] The self-assembly of G x M y H z BCPs is driven not only by the hydrophilic-hydrophobic balance of the block copolymers but also by the non-covalent interactions (pi-pi and H-bonding) between the GFF sequences, more specifically due to the FF moieties. As the self-assembly based on FF is stabilized mostly via H bonding (FF-FF and FFsolvent) which is weakened or even completely suppressed at high temperature, [START_REF] Wang | Peptide self-assembly: thermodynamics and kinetics[END_REF][START_REF] Huang | Temperatureinduced reversible self-assembly of diphenylalanine peptide and the structural transition from organogel to crystalline nanowires[END_REF] the morphologies observed at ambient temperature (about 25-30 °C) might not be identical with the morphologies formed during the polymerization at 70 °C. In our previous work, the G 65 M 7 H 28 suspension was shown to be composed of worm-like structures at 70 °C, and of flower-like fiber bundles at 20 °C. 1 Therefore, the effect of temperature on the self-assembly of the G 65 M 7 H z series was investigated. As expected, at 70 °C, no fibers or fiber bundles were observed for any G 65 M 7 H z BCP. Instead, pure worms, a mixture of spheres and worms, beaded worms ( pearl necklace-like morphology) and a mixture of worms and spheres were observed for G 65 M 7 H 28 , G 65 M 7 H 42 , G 65 M 7 H 54 , and G 65 M 7 H 63 , respectively (Fig. 5). In contrast, at both 70 °C and 30 °C, G 65 M 7 H 200 formed spheres (Fig. 5). After a critical block length, high temperature has no effect on the self-assembly since the impact of the SAP becomes negligible.

This temperature-induced morphological change was also detected by DLS (Fig. 6). The suspension of G 65 M 7 H 42 at 20 °C presented a broad average hydrodynamic diameter distribution, in agreement with the highly anisotropic structure observed in TEM. Upon heating from 20 °C to 70 °C, the average hydrodynamic particle diameter distributions gradually shifted towards smaller dimensions, and above 60 °C, only submicron particles were detected. This change was completely reversible as shown by the results of the DLS measurements carried out at 20 °C before and after heating to 70 °C, which were essentially identical. In contrast, the DLS data of G 65 M 7 H 200 were not affected by heating up to 70 °C (Fig. 6). This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

The particle hydrodynamic diameter distributions remained centred around 100 nm and became slightly narrower upon heating from 20 to 70 °C (Fig. 5), in agreement with the TEM images.

Although heating disrupts the associative interactions between the GFF moieties, low temperatures have been reported to increase the hydrophilicity of PHPMA. [START_REF] Lovett | Order-Order Morphological Transitions for Dual Stimulus Responsive Diblock Copolymer Vesicles[END_REF] Previous studies demonstrated that PGMA-b-PHPMA BCP morphologies undergo a worm-to-sphere transition upon cooling from 20 °C to 5 °C. [START_REF] Verber | Rheological studies of thermo-responsive diblock copolymer worm gels[END_REF][START_REF] Blanazs | Sterilizable Gels from Thermoresponsive Block Copolymer Worms[END_REF] Here, cooling also affected the P G x M y H z BCP morphologies. Upon cooling from 30 °C to 4 °C the self-assembled structures of the G 65 M 7 H z series changed markedly from defined bundles to spheres. In the case of G 65 M 7 H 28 , cooling transformed the fiber bundles into not very defined fluffy bundles/balls (G 65 M 7 H 28 ). In the case of G 65 M 7 H 42 , a temperature decrease led to the formation of individual spheres and small bundles of fibers, while for the G 65 M 7 H 54 case the temperature change resulted in the formation of spheres coated with short fibers (sea urchins). For the longest hydrophobic block (G 65 M 7 H 200 ) there was no transformation as the observed morphology at 30 °C was spheres. These unique structures have never been reported before once again highlighting the powerful structuring ability of the SAP moieties.

P(GMA-stat-(MAm-GFF))-b-PHPMA via RAFT dispersion PISA in water/ethanol mixtures

As mentioned above, some of the G x M y H z suspensions were not colloidally stable, likely due to the interactions between the GFF moieties located in the periphery of the selfassembled structures. To deal with this instability, RAFT polymerization was conducted in a water/ethanol binary mixture instead of pure water (MAm-GFF does not dissolve in water while it dissolves in ethanol up to 5% w/w at 30 °C). The RAFT dispersion polymerisation of HPMA using mCTA 2 (G 65 M 7 ) was carried out in 50/50, 60/40 and 80/20 water/ ethanol mixtures. These polymerization reactions reached 100% conversion of HPMA (G 65 M 7 H 66 ), forming stable colloidal suspensions over the course of 12 months. Although the three reaction mixtures were slightly opaque, the samples in ethanol-rich binary mixtures (50/50 and 60/40) did not form defined structures as ethanol is a good solvent for PHPMA. In contrast, in an 80/20 water/ethanol mixture short rigid fibres of about 500 nm length were observed (Fig. 7). Unlike the fibrils formed in pure water (Fig. 4), these short fibers did not bundle.

Conclusions

This work focuses on the morphology of P(GMA-stat-(MAm-GFF))-b-PHPMA (GMH) diblock copolymers synthesized via a one-step RAFT-mediated PISA approach. The effects of the composition of the GFF-containing corona-forming block and of the length of the core-forming block on the colloidal stability and the morphologies of the self-assembled structures were examined. Some of the PISA morphologies were colloidally stable for months while others sedimented and phase separated within a few hours after or even during the polymerization. Block extension of the shortest mCTA with the highest ratio of MAm-GFF/GMA units (mCTA 1) resulted in the formation of bundled fibres, similar to the morphology of mCTA1 alone. When a longer mCTA with a lower percentage of MAm-GFF (mCTA 2) was used, bundled fibers were again observed for low DPs of PHPMA. However, for a DP of PHPMA of above 54, first urchin-like structures and ultimately spheres were formed. This transformation suggests that above a threshold DP of the core-forming PHPMA block, the effect of the SAP becomes negligible and the self-assembly follows the "classical" regime dictated by the hydrophilic/hydrophobic balance under PISA conditions. Further reducing the MAm-GFF content in the mCTA (mCTA 3) resulted in the formation of urchin-like morphologies and finally spheres when the DP of PHPMA reached 510. All the examined formulations resulted in the formation of defined morphologies at 70 °C, 30 °C and 4 °C. At 70 °C, the impact of the SAP was weakened as only worms and spheres were observed (H-bonds are disrupted at high T ). Upon cooling to ambient temperature, bundles of fibers were found as the main structures. The formation of such bundled fibres is the fingerprint of the GFF-sequence. Upon further cooling down to 4 °C, the morphology transformed from fiber bundles to spheres coated with short fibers that resemble sea urchins. Finally, the PISA of G 65 M 7 H 66 in an 80/20 water/ethanol binary mixture resulted in the formation of colloidally stable non-bundled short fibers. This study confirms that including SAP moieties within the structure of block copolymers is an elegant approach to prepare self-assembled structures with rare and complex morphologies. Such a strategy provides a powerful tool for structural derivatization that can be used in combination with different SAP sequences and polymerization techniques. The reported novel morphologies could find use in the ever-changing fields of dynamic materials and drug delivery.

Fig. 1

 1 Fig.1Representative TEM images of mCTA 1 (G 24 M 3 ) and mCTA 2 (G 65 M 7 ) in Milli-Q water at 0.1 wt% and prepared at ambient temperature (about 25 °C). Note: no morphology was detected for mCTA 3.

Fig. 3

 3 Fig. 3 Photographs of the G 24 M 3 H 14 (left) and G 24 M 3 H 58 (right) BCP reaction media after 24 h (aqueous RAFT-PISA carried out at 10 wt% in water at 70 °C).

Fig. 2

 2 Fig.2DMF SEC data (refractive index detector) for the P(GMA 24 -stat-(MAm-GFF) 3 ) mCTA (mCTA 1) and the diblock copolymers G 24 M 3 H 14-108 .

Fig. 4

 4 Fig. 4 Representative TEM images of G 24 M 3 H 14-108 (top row), G 65 M 7 H 42-200 (middle row) and G 200 M 9 H 51-510 (bottom row) diblock copolymers prepared by RAFT-PISA in water at 10 wt% solids. Samples were diluted 100-fold at 30 °C. Red contour indicates the colloidally stable samples.

Fig. 5

 5 Fig. 5 Representative TEM images of G 65 M 7 H 28-200 diblock copolymers at 70 °C (top row), 30 °C (middle row) and 4 °C (bottom row) prepared by RAFT-PISA in water at 10 wt% solids. The samples were prepared at the corresponding temperatures (100-fold dilution and deposition on the TEM grid).

Fig. 6

 6 Fig. 6 Evolution with temperature of the hydrodynamic diameter of G 65 M 7 H 42 (top) and G 65 M 7 H 200 (bottom) BCP objects at 0.1 wt% in Milli-Q water.

Table 1

 1 Molecular characterization of P(GMA-stat-(MAm-GFF)) macro-CTAs

	Name	Abbreviation a	M n NMR (g mol -1 )	M n GPC (DMF) (g mol -1 )	Đ GPC (DMF)
	mCTA 1	G 24 M 3	5160	6600	1.35
	mCTA 2	G 65 M 7	13 470	13 300	1.28
	mCTA 3	G 200 M 9	35 970	32 100	1.36

Scheme 1 Schematic representation of the synthesis of the SAP-containing methacrylamide (MAm-GFF) and P(GMA x -stat-(MAm-GFF) y )-b-PHPMA z diblock copolymers.PaperPolymer Chemistry

Acknowledgements

The PDRA fellowship of T. P. Tuyen Dao was managed by the French National Research Agency (ANR): "investissement d'Avenir" LabEX CheMISyst, grant number ANR-10-LABX-05-01. The authors also thank CNRS for funding this work via the "Osez l'interdisciplinarité-2017" programme awarded to M. Semsarilar.

Conflicts of interest

There are no conflicts to declare. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence.

Notes and references