Stéphanie Chevalier

Vincent Noël

Laurence Calzone

Andrei Zinovyev

Loïc Paulevé

Synthesis and Simulation of Ensembles of Boolean Networks for Cell Fate Decision

The construction of models of biological networks from prior knowledge and experimental data often leads to a multitude of candidate models. Devising a single model from them can require arbitrary choices, which may lead to strong biases in subsequent predictions. We introduce here a methodology for a) synthesizing Boolean model ensembles satisfying a set of biologically relevant constraints and b) reasoning on the dynamics of the ensembles of models. The synthesis is performed using Answer-Set Programming, extending prior work to account for solution diversity and universal constraints on reachable fixed points, enabling an accurate specification of desired dynamics. The sampled models are then simulated and the results are aggregated through averaging or can be analyzed as a multi-dimensional distribution. We illustrate our approach on a previously published Boolean model of a molecular network regulating the cell fate decisions in cancer progression. It appears that the ensemble-based approach to Boolean modelling brings new insights on the variability of synergistic interacting mutations effect concerning propensity of a cancer cell to metastasize.

Introduction

The ability to derive one single model from observations of a biological system usually faces arbitrary choices, sometimes referred to as art.

Computational models of molecular interaction networks are usually built from data related to the architecture of the network from known interactions; and data related to its dynamics, such as measurements of gene expressions or proteins activity at different times and/or conditions. However, despite huge advances in experimental technologies, observations of the biological processes stay very scarce, either in terms of temporal resolution, number of observed entities, synchronisation between measure points, or a variety of experimental conditions. Combined with complex structures for molecular interactions, the model engineering problem, in this case, appears to be largely under-specified, leading to (too) many potential candidate models.

co-first authors

Boolean Networks (BNs), and logical models in general, are widely adopted for the modelling of signalling pathways and gene and transcription factors networks [START_REF] Collombet | Logical modeling of lymphoid and myeloid cell specification and transdifferentiation[END_REF][START_REF] Zañudo | Discrete dynamic network modeling of oncogenic signaling: Mechanistic insights for personalized treatment of cancer[END_REF][START_REF] Clarke | Executable cancer models: successes and challenges[END_REF]. With BNs, the activity of components is caricatured to "off" and "on", and their evolution is computed according to logical rules (e.g., gene 1 can be active only whenever its activators 2 and 3 are active). However, in practice, biological data still let open a multitude of candidate BNs. Thus, arbitrary modelling choices have to be made, e.g., by prioritizing certain logics between regulators or by preferring smallest/largest models, which may introduce biases in subsequent model predictions.

In this paper, we present an approach aiming at reducing modelling biases by constructing and reasoning on dynamics of ensembles of BNs. The idea of ensemble modelling has recently gained momentum with machine learning, notably with random forests. By analogy, we constitute ensembles of BNs sampled from the whole multitude of models compatible with network architecture and dynamical properties. They are then simulated asynchronously and the simulations are aggregated through averaging. The obtained results allow an interpretation at the level of cell population and take into account its potential heterogeneity.

In the literature, ensembles of random BNs have been employed to show emerging properties of families of BNs sharing properties related to their architecture or logic rules [START_REF] Kauffman | A proposal for using the ensemble approach to understand genetic regulatory networks[END_REF][START_REF] Krawitz | Basin entropy in boolean network ensembles[END_REF]. In [START_REF] Schwieger | Graph representations of monotonic boolean model pools[END_REF], ensembles of BNs sharing a network architecture are used to assess dynamical properties of qualitative differential equations. In contrast, our approach is focused on ensembles of models which satisfy a set of constraints both on their architecture and on their dynamics. Synthesis of BNs from such constraints received a lot of interest in the literature [START_REF] Corblin | A declarative constraint-based method for analyzing discrete genetic regulatory networks[END_REF][START_REF] Terfve | CellNOptR: a flexible toolkit to train protein signaling networks to data using multiple logic formalisms[END_REF][START_REF] Dorier | Boolean regulatory network reconstruction using literature based knowledge with a genetic algorithm optimization method[END_REF][START_REF] Yordanov | A method to identify and analyze biological programs through automated reasoning[END_REF][START_REF] Goldfeder | BRE:IN -a backend for reasoning about interaction networks with temporal logic[END_REF][START_REF] Chevalier | Synthesis of Boolean Networks from Biological Dynamical Constraints using Answer-Set Programming[END_REF], and methods like [START_REF] Corblin | A declarative constraint-based method for analyzing discrete genetic regulatory networks[END_REF][START_REF] Yordanov | A method to identify and analyze biological programs through automated reasoning[END_REF][START_REF] Goldfeder | BRE:IN -a backend for reasoning about interaction networks with temporal logic[END_REF][START_REF] Chevalier | Synthesis of Boolean Networks from Biological Dynamical Constraints using Answer-Set Programming[END_REF] allows reasoning implicitly on ensembles of models, notably by enabling checking for their emptiness.

Here, we extend prior work on BNs synthesis from reachability and attractor properties with Most Permissive semantics [START_REF] Chevalier | Synthesis of Boolean Networks from Biological Dynamical Constraints using Answer-Set Programming[END_REF] to support universal properties on (reachable) fixed points and the specification of network perturbations. The synthesis is performed with the logic framework of Answer-Set Programming [START_REF] Baral | Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF] (ASP). Then, we use heuristics to drive the ASP solver in different regions of the solution space to sample ensembles of BNs capturing the diversity of the comprehensive solution set. Dynamics of ensembles are then explored through stochastic simulations for quantifying the propensities of reachable attractors, subject to different network perturbations. To that aim, we extended the simulator MaBoSS [START_REF] Stoll | MaBoSS 2.0: an environment for stochastic boolean modeling[END_REF] to support ensembles of BNs as input.

We illustrate our approach on a model of molecular pathways regulating tumour invasion and migration [START_REF] Cohen | Mathematical modelling of molecular pathways enabling tumour cell invasion and migration[END_REF]. We sampled ensembles of BNs sharing the same network architecture as the original model and constrained by the dynamical properties related to attractor reachability. Then, as in the original study, we evaluated the shift of reachable phenotypes caused by an epistatic interaction between mutations in model genes (gain of Notch function and loss of p53 function). It appears that, contrary to the initial single model analysis, the ensemble approach reveals a potential variability in the effectiveness of the double mutant to enhance the metastasis potential.

Background

Boolean Networks

A Boolean network (BN) of dimension n is a function

f : B n → B n (1)
where B := {0, 1}. For all i ∈ [n], f i : B n → B denotes the local function of the i-th component. A vector x ∈ B n is called a configuration of the BN f . The set of components which value differs between two configurations x, y ∈ B n is denoted by ∆(x, y)

:= {i ∈ [n] | x i = y i }.
A BN f is said locally monotonic whenever each of its local functions is monotonic (this does not imply f monotonicity). Intuitively, local monotonicity imposes that a variable always appears with the same sign in a minimal disjunctive/conjunctive normal form of the local functions.

Fig. 1 is an example of locally-monotonic BN with n = 3.

Mutations

In the following, we will consider the analysis of a BN f subject to some permanent perturbations of its components, that we refer to as mutations, being either a gain of function (GoF; locked to 1) and loss of function (LoF; locked to 0). A mutation is specified by a couple (i, v), where i ∈ [n] is a component and v ∈ B is its forced value. Given a BN f and a set of mutations M ⊆ [n] × B, we denote by f /M the mutated BN, where, for each i

∈ [n], (f /M ) i (x) := v if (i, v) ∈ M , and (f /M ) i (x) := f i (x) otherwise.
Influence graph For each i ∈ [n], f i typically depends on a small subset of components of the BN. The influence graph summarizes these dependencies with a positive (resp. negative) edge from node j to i if there are configurations in which the sole increase of j would strictly increase (resp. decrease) the value of f i . A node can have both positive and negative influences on i, indicating that f i is non-monotonic. Remark that different BNs can have the same influence graph. Fig. 1 (right) shows the influence graph of the BN example.

Definition 1. Given a BN f of dimension n, its influence graph G(f ) is a directed graph ([n], E + , E -) with positive and negative edges such that (j, i) ∈ E + (resp. (j, i) ∈ E -) iff ∃x, y ∈ B n s.t. ∆(x, y) = {j}, x j < y j , and

f i (x) < f i (y) (resp. f i (x) > f i (y)). The influence graph G = ([n], E + , E -) is a subgraph of G = ([n], E + , E -), denoted by G ⊆ G , iff E + ⊆ E + and E -⊆ E -. f1(x) := ¬x2 f2(x) := ¬x1 f3(x) := ¬x1 ∧ x2

BN Semantics

From a configuration x ∈ B n , semantics of BNs specify how to compute the next possible configurations. One of the most classical semantics is the fullyasynchronous (often simply called asynchronous), where only one component i is updated at a time (to the value f i (x)). It can be defined as a binary relation However, as demonstrated in [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF], the fully-asynchronous semantics of BNs, as the synchronous and (general) asynchronous, are not faithful abstractions of quantitative systems: they can both introduce spurious behaviours (as expected with qualitative models) and miss others.

∀x, y ∈ B n , x f -→ a1 y iff ∃i ∈ [n] : ∆(x, y) = {i} ∧ y i = f i (x) .
The Most Permissive (MP) semantics of BNs [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF] offers the guarantees to not preclude any behaviour realisable in any quantitative refinement of the model, thus providing a formal over-approximation of dynamics. Moreover, the abstraction is minimal: any behaviour it predicts is realisable by a quantitative refinement of the BN using the asynchronous semantics. Importantly, the complexity for deciding main dynamical properties is considerably lower than with (a)synchronous semantics, as we will mention in the next subsection.

MPBNs can be defined by the means of hypercubes (partially) closed by f , a hypercube being specified by a vector associating each component to either a fixed Boolean value or free ( * ).

Definition 3 (Hypercube).

A hypercube h of dimension n is a vector in (B ∪ { * }) n . The set of its associated configurations is denoted by c(h

) := {x ∈ B n | ∀i ∈ [n], h i = * ⇒ x i = h i }. Given two hypercubes h, h ∈ (B ∪ { * }) n , h is smaller than h if and only if ∀i ∈ [n], h i = * ⇒ h i = h i . Definition 4 (K-closed hypercube). Given a subset of components K ⊆ [n], a hypercube h ∈ (B ∪ { * }) n is K-closed by f whenever for each configuration x ∈ c(h), for each component i ∈ K, h i ∈ { * , f i (x)}. It is minimal whenever no different K-closed hypercube is smaller than it.
A hypercube [n]-closed by f is also known as a trap space [START_REF] Klarner | Computing maximal and minimal trap spaces of boolean networks[END_REF]. Starting from a configuration x ∈ B n , the MP semantics allows transitions towards any configuration y which is present in at least one smallest K-closed hypercube h containing x, for some K ⊆ n, and so that the state of each component i ∈ K of y can be computed by f i from a configuration of h.

Example 1. Let us consider the BN f : B 3 → B 3 with f 1 (x) := 1, f 2 (x) := x 1 , and f 3 (x) := x 1 ∧ ¬x 3 . The hypercube 1 * * is closed by f , with c(1 * * ) = {100,
Definition 5 (Most-Permissive semantics). Given a BN f of dimension n and two configurations x, y ∈ B n , y ∈ ρ f mp (x) if and only if there exists K ⊆ [n] such that the smallest K-closed hypercube h containing x verifies (1) y ∈ c(h), and (2) ∀i ∈ K, there exists a configuration z ∈ c(h) such that f i (z) = y i .
Fig. 2 gives examples of computations of ρ f mp . A way to interpret the MP semantics is to see the components free in a hypercube as being in the course of changing of state, while other components can independently consider them either as 0 or 1. This abstracts the missing information on the ordering of thresholds of activation/inhibition between components: while the quantitative value of component u progressively increases, at a given time it can be high enough to activate a component (i.e., [START_REF] Baral | Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF] but not yet high enough to activate another one (i.e., 0). These dynamic states are overlooked by asynchronous semantics, making it an incorrect over-approximation of quantitative systems, contrary to the MP semantics [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF].

Dynamical properties

In the following, we will focus on two main dynamical properties of BNs: reachability which relates to the existence of trajectories between two configurations, and attractors which relates to long-run behaviours by identifying the smallest sets of configurations closed by reachability. Definition 6 (Reachability). Given two configurations x, y ∈ B n of a BN f with semantics σ, y is reachable from x whenever y ∈ ρ f σ (x).

Definition 7 (Attractor). A non-empty set of configurations A ⊆ B n is an attractor of the BN f with semantics σ whenever ∀x ∈ A, ρ f σ (x) = A. When A = {x} for some x ∈ B n , we say that x is a fixed point.

With MP semantics, attractors match with minimal trap spaces. With fullyasynchronous semantics, deciding if y ∈ ρ f a1 (x) or if x belongs to an attractor are both PSPACE-complete problems. With MP semantics, deciding if y ∈ ρ f mp (x) is PTIME if f is locally-monotonic and P NP otherwise; deciding if x belongs to an attractor is coNP if f is locally-monotonic and coNP coNP otherwise. Deciding if there exists a fixed point is NP-complete with both semantics [START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF].

Notice the following relations between MP and fully-asynchronous semantics:

x ∈ B n is a fixed point with MP semantics if and only if it is a fixed point with full-asynchronous semantics (iff f (x) = x); y ∈ B n is reachable from x ∈ B n with the fully-asynchronous semantics only if it is reachable with MP semantics (ρ f a1 (x) ⊆ ρ f mp (x)); the number of attractors with MP semantics is less than or equal to the number of attractors with fully-asynchronous semantics.

Answer-Set Programming

Answer Set Programming (ASP; [START_REF] Baral | Knowledge Representation, Reasoning and Declarative Problem Solving[END_REF][START_REF] Gebser | Answer Set Solving in Practice[END_REF]) is a declarative approach to solving combinatorial satisfaction problems. It is close to SAT (propositional satisfiability) [START_REF] Lin | ASSAT: Computing answer sets of a logic program by SAT solvers[END_REF] and known to be efficient for enumerating solutions of NP problems comprising up to tens of millions of variables while providing a convenient language for specifying the problem. We give a very brief overview of ASP syntax and semantics that we use in the next sections; see [START_REF] Gebser | Answer Set Solving in Practice[END_REF] for more details. An ASP program is a Logic Program (LP) being a set of logical rules with first order logic predicates of the form:

1 a0 ← a1, . . ., an, not an+1, . . ., not a n+k .

where a i are (variable-free) atoms, i.e., elements of the Herbrand base, which is built from all the possible predicates of the LP. The Herbrand base is built by instantiating the LP predicates with the LP terms (constants or elements of the Herbrand universe).

Essentially, such a logical rule states that when all a 1 , . . . , a n are true and none of a n+1 , . . . , a n+k can be proven to be true, then a 0 has to be true as well. Whenever a 0 is ⊥ (false), the rule, also called integrity constraint, becomes: 2 ←a1, . . ., an, not an+1, . . ., not a n+k . Such a rule is satisfied only if the right-hand side of the rule is false (at least one of a 1 , . . . , a n is false or at least one of a n+1 , . . . , a n+k is true). On the other hand, a0 ← (a 0 is always true) is abbreviated as a0.

A solution (answer set) is a stable Herbrand model, that is, a minimal set of true atoms where all the logical rules are satisfied. For instance, consider the following program:

a. b ← a. d ← a, c.
It has for unique solution {a, b}: indeed, whereas {a, b, d} does not contradict the rules, d is not a fact and cannot be derived from a rule; so it is not stable.

ASP allows using variables (starting with an upper-case) instead of terms/predicates: these template declarations will be expanded before the solving. We also use the notation a(X): b(X) which is satisfied when for each b(X) true, a(X) is true. If any term follows such a condition, it is separated with ;.

ASP can express disjunctive logic programs [START_REF] Lobo | Foundations of disjunctive logic programming[END_REF], by the means of disjunctions in the rule head (";"-separated atoms):

a; b ← body.
Such a disjunctive rule implies that solutions are subset minimal: an answer set is a solution only if none of its subsets is itself a solution [START_REF] Eiter | Answer Set Programming: A Primer[END_REF]. For instance, let's consider the disjunction: a; b; c.

The interpretation I = {a, b} is a model but not minimal: both interpretations {a} and {b} are smaller than I and satisfy the rule. Hence I is not a solution. As showed in [START_REF] Eiter | On the computational cost of disjunctive logic programming: Propositional case[END_REF], the complexity of problems addressed with ASP can be extended thanks to disjunctive rules up to 2QBF, i.e. a two quantification levels Boolean formula (∀x∃y.φ or ∃y∀x.φ where φ is a quantifier-free propositional formula). Indeed, 2QBF can be reduced to the problem of verifying the existence of an answer set of a disjunctive ASP program.

BN Synthesis from Architecture and Dynamical Properties

We formulate the problem of BN synthesis as a Boolean satisfiability problem encoded in ASP. With this approach, we leverage a priori knowledge and experimental data as constraints on the network architecture and the dynamical properties of the models under the MP semantics. Our method is based on [START_REF] Chevalier | Synthesis of Boolean Networks from Biological Dynamical Constraints using Answer-Set Programming[END_REF], which implements constraints on existence and absence of trajectories between partially-specified configurations, existence of (reachable) fixed points and trap spaces. In biological applications, these constraints match well the observed properties of cell populations evolving towards mutually exclusive phenotypes.

In this paper, we extend [START_REF] Chevalier | Synthesis of Boolean Networks from Biological Dynamical Constraints using Answer-Set Programming[END_REF] to support universal properties on (reachable) attractors. This enables specifying tight dynamical constraints. For instance, given a set of experimentally observed phenotypes, existential constraints guarantee that at least one attractor of the model dynamics match with each phenotype, whereas a universal constraint ensures that every attractor matches with at least one of the phenotype.

A universal property involves by nature universal quantifiers. ASP can address formulas implying one level of universal quantifier (i.e., of the form ∃x∀y : P (x, y)) thanks to a technique presented in [START_REF] Eiter | On the computational cost of disjunctive logic programming: Propositional case[END_REF]. To explore a set of values and check the respect of a property for each, it uses a disjunctive rule and a saturation on the same term. A disjunctive rule implies the subset minimality semantics. This minimality ensures an answer set is a solution only if none of its subsets is itself a solution [START_REF] Eiter | Answer Set Programming: A Primer[END_REF]. Hence, saturating the answer set with the predicates of the disjunction cleverly exploits this minimality: the solver is forced to ensure that no strict subset of these predicates form a solution.

Universal constraints on fixed points

We exploit this saturation technique [START_REF] Eiter | On the computational cost of disjunctive logic programming: Propositional case[END_REF][START_REF] Eiter | Answer Set Programming: A Primer[END_REF] for ensuring universal constraints on the fixed points or fixed points reachable from a given configuration. We describe here the ASP rules for the universal fixed point constraint, which ensures that all the fixed points of the BN are compatible with a given set of markers (observations). To that purpose, we let the solver deduce a configuration z by the disjunctive rule: cfg(z,N,-1) ; cfg(z,N,1) ← node(N).

The predicate template cfg(X,N,V) assigns the value V to the literal N in the configuration X. Through the above rule, a set of node values is thus constituted to define a configuration z, with the predicate cfg(z,N,_) subject to the subset minimality semantics. To respect the desired property, each configuration z is either not a fixed point (f (z) = z) or has the same component states than the ones expressed in a dedicated predicate. A configuration is not a fixed point whenever at least one of its component can change of state: mcfg(z,N,V) ← cfg(z,N,V). valid ← cfg(z,N,V) ; eval(z,N,-V). mcfg(X,N,V) predicate template leads to the evaluation of the configuration X given the Boolean rules of the network [START_REF] Chevalier | Synthesis of Boolean Networks from Biological Dynamical Constraints using Answer-Set Programming[END_REF]. The reachable values are then stored in the predicate eval(X,N,V). Whenever it is possible to evaluate a component N to the opposite value than in z, then z is not a fixed point, making valid true.

Otherwise, z has to have the same component states than those specified by an observation X marked by the predicate is_universal_fp(X), which is expressed by the following ASP rule: valid ← cfg(z,N,V):obs(X,N,V); is_universal_fp(X).

Observe in l.10 and l.11 that each time an assignment is in agreement with the desired property, a predicate valid is deduced, which triggers the saturation of the configuration z:

cfg(z,N,-V) ← cfg(z,N,V), valid.
Thus, when valid is deduced, the answer set contains all possible component values for z. According to the subset minimality semantics, the solver is then forced to ensure that there is no sub-answer set. And the only way to find such a smaller answer-set is to find a z from which valid cannot be deduced, i.e., which is a counter-example to the universal property: in that case, l.13 eliminates the answer set:

13 ← not valid.

A variant of this constraint enables to restrict the universal property to fixed points that are reachable from a given initial configuration. This is specified by is_universal_fp(X,S) predicates, where S points the initial configuration, and X to an observation, as used above. By combining such predicates, one can then specify sets of phenotypes reachable from a given configuration. The encoding of this variant contains a third way to deduce valid: the non-reachability of the configuration z from S.

Our implementation also offers to specify mutations, which can be combined with reachability and with universal constraints on reachable fixed points to leverage observations about cell fates in different mutation conditions.

Synthesis problem

Synthesis requires (i) an influence graph to delimit the interactions that can be used by the BNs and (ii) the dynamical properties of the behaviours that have to be reproduced. For modelling the tumour invasion and migration as in [START_REF] Cohen | Mathematical modelling of molecular pathways enabling tumour cell invasion and migration[END_REF], the dynamical properties refer to cell fate observations in different mutation conditions. These fates are described by sets of markers (i.e. a set of values for some nodes of the network) which constitutes partial observations of genes activity. In term of dynamics, these observations are related to reachable attractors in the corresponding mutated BNs.

A (partial) observation o of a configuration of dimension n is specified by a set of couples associating a component to a Boolean value: o ⊆ [1]n×B, assuming there is no i ∈ [n] such that {(i, 0), (i, 1)} ⊆ o.

Formally, the synthesis problem we tackle is the following. Given find a locally-monotonic BN f of dimension n such that -G(f ) ⊆ G, there exist p configurations x 1 , . . . , x p such that:

• (observations) ∀m ∈ [p], ∀(i, v) ∈ o m , x m i = v, • (positive reachability) ∀(m, m ) ∈ PR, x m ∈ ρ f mp (x m ), • (fixed points) ∀m ∈ FP, f (x m ) = x m , • (universal fixed point) ∀z ∈ B n , f (z) = z ⇒ ∃m ∈ UFP : ∀(i, v) ∈ o m , z i = v; • (universal reachable fixed point) ∀(q, m ) ∈ URFP, ∀z ∈ ρ f mp (x q ), z = f (z) ⇒ ∃(q, m) ∈ URFP : ∀(i, v) ∈ o m , z i = v.
Each of these constraints can be parametrized by mutations, in which case, the properties have to be verified on the mutated f . Remark that such a problem can be non-satisfiable.

Our encoding also offers constraints related to the absence of paths between configurations (negative reachability) and to trap space where a set of components have a fixed state matching with a given observation [START_REF] Chevalier | Synthesis of Boolean Networks from Biological Dynamical Constraints using Answer-Set Programming[END_REF]. Moreover, one can optionally impose that the influence graph of f is equal to the input G.

Our implementation avoids redundancy in the models by enumerating only among non-equivalent BNs (i.e., their values differ for at least one configuration). This is achieved by using a canonical representation of Boolean functions in disjunctive normal form with a total ordering between clauses.

In total, our encoding generates O(ndk 2 ) atoms and O(nd 2 k 2 ) rules, where d is the in-degree of nodes in the influence graph, and k is a fixed bound on the number of clauses of Boolean functions. Whenever k is set to d d/2 , the complete set of solutions can be enumerated.

Sampling the diversity of all solutions

The whole set of constraints, comprising the domain of admissible BNs and the dynamical properties they should satisfy, is represented by a single logic program expressed in ASP, such that each solution corresponds to a distinct BN.

Whereas the enumeration of ASP solutions is known to be efficient, typical solvers will enumerate solutions by slightly varying parts of a firstly identified one. Thus, a partial enumeration will very likely give a set of solutions which are all look alike, e.g., where the Boolean function of only one component varies.

Inspired by [START_REF] Razzaq | Computing diverse boolean networks from phosphoproteomic time series data[END_REF], we tweak heuristics of the solver clingo [START_REF] Gebser | Clingo = ASP + control: Preliminary report[END_REF] to stir it towards distant solutions: at each solution, we randomly select a subset of variables assignments and ask the solver to avoid them in the next iterations. At the cost of enumeration speed, this allows sampling ensembles of diverse BNs.

Stochastic Simulations of Ensembles of BNs

Continuous-time Boolean modelling

We first recall the continuous-time Markov chain interpretation of BNs introduced in [START_REF] Stoll | Continuous time boolean modeling for biological signaling: application of gillespie algorithm[END_REF]. Considering a BN f of dimension of n, we represent the state evolution by a Markov process s : t → s(t) defined on t ∈ I ⊂ R applied on the network state space, with I the simulation interval. This process is defined by: 1. Its initial condition:

P [s(0) = x], ∀x ∈ B n 2.
Its conditional probabilities (of a single condition):

P [s(t) = y|s(t ) = x], ∀x, y ∈ B n , ∀t , t ∈ I, t < t
In continuous-time, these conditional probabilities are defined as transition rates [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF]: ρ(x → y)(t) ∈ [0, ∞]. Because we want a generalization of the fullyasynchronous Boolean dynamics, transition rates ρ(x → y) are non-zero only

if x f -→ a1 y, i.e., a single component i ∈ [n] is changing of value. In that case, each local function f i (x) is replaced by two functions R up/down i (x) ∈ [0, ∞].
The transition rates are defined as follows:

ρ(x → y) = R up i (x) if x i = 0 R down i (x) if x i = 1 with ∆(x, y) = {i}
where R up i corresponds to the activation rate of node i, and R down i corresponds to the inactivation rate of node i. Therefore, the continuous Markov process is completely defined by all these R up/down and an initial condition. By default, the value of these rates is set to 1, but they can be modified to represent the time scales of different processes.

To explore the probability space of this Markov process, we use the Gillespie algorithm [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF]. This algorithm produces a set of realizations or stochastic trajectories of the Markov process. From this finite set, probabilities can be estimated.

To relate continuous-time probabilities to real processes, an observable time window δt is defined. A discrete-time τ ∈ N stochastic process s(τ ) can be extracted from the continuous-time Markov process:

P [s(τ ) = x] ≡ 1 δt (τ +1)δt τ δt P [s(t) = x]dt
For each trajectory j, we compute the time for which the system is in state x in the window [τ δt, (τ + 1)δt], and divide it by δt. We obtain an estimate of P [s(τ ) = x] for trajectory j, i.e. Pj [s(τ ) = x]. Then to compute the estimate of a set of trajectories, we compute the average over j of all Pj [s(τ ) = x].

Lifting to ensembles of BNs

To simulate an ensemble of BNs, we first choose a total number of stochastic trajectories M . We generate M/k stochastic trajectories for each model and compute the average Pk [s(τ ) = x] for all models k. We then compute the average over k of all Pk [s(τ ) = x], to obtain the P [s(τ ) = x] for the ensemble of boolean networks. We also keep the option to export the individual probability distributions Pk [s(τ ) = x] to allow us analyzing the composition of the ensemble. The approach results in time-series of the probability for each observed state. The case study hereafter focuses on steady-state analysis. This imposes to simulate the ensemble long enough to reach stationarity, requiring a preliminary analysis. We can then study the proportion of each attractor for our ensemble.

We implemented this new feature in the MaBOSS simulation software [START_REF] Stoll | MaBoSS 2.0: an environment for stochastic boolean modeling[END_REF]

Case Study on Cell Fate Decision Modelling

Background Model

We illustrate our ensemble modelling approach on a published model of cell fate decision leading to the early events of the metastasis or cell death through apoptosis [START_REF] Cohen | Mathematical modelling of molecular pathways enabling tumour cell invasion and migration[END_REF]. Initial triggers, such as DNA damage or micro-environmental cues, and the activity of some genes or proteins participating in the process affect the final decision. The signalling pathway involves TGFbeta, WNT, beta-catenin, p53 and its homologs, selected miRNA, and transcription factors of the epithelial to mesenchymal (EMT) transitions. Fig. 3 shows the influence graph of the BN.

The functions of the BN, we refer to as "Cohen's model", have been designed manually so the simulations fit with experimental data related to stable phenotypes under different single mutations. Then, the initial publication explored the synergy between mutants that led to metastatic phenotypes.

Single Model Analysis

We first reproduced part of the analysis of [START_REF] Cohen | Mathematical modelling of molecular pathways enabling tumour cell invasion and migration[END_REF] on the original Cohen's model by computing the propensities of attractors reachable from 4 possible initial conditions, where all nodes are inactive, except miRNAs that are active, and the 2 nodes modelling DNA damage and micro-environment cues that are free. We considered the wild-type condition (Fig. 4(a)) with no mutation, and the double-mutant of p53 LoF and Notch GoF (Fig. 4(b)).

The wild-type model has 9 fixed points that each correspond to one of the 4 identified physiological phenotypes: Apoptosis, EMT, Metastasis (or equivalent to Migration) and Homeostatic State (HS). The double-mutant shows exclusively the Metastasis phenotype.

Ensemble Analysis

Synthesis To test the impact of alternative Boolean functions, we synthesised ensembles of BNs that share the same influence graph as Cohen's model and reproduce the desired dynamics. We synthesized two ensembles of 1,000 diverse BNs each, where we disallowed having cyclic attractors. The first ensemble ensures only the wild-type (WT) behaviour, meaning that all the fixed points match with one of the 4 physiological phenotypes, and each physiological phenotype is reachable from at least one of the initial condition. The second ensemble adds further constraints related to the single mutations of p53 LoF, which should show the same behaviour as WT, and Notch GoF, where only 2 of the WT phenotypes and a third different one should be observed 6 .

Ensemble Simulations With the same settings as with Cohen's model, we performed stochastic simulations of the two synthesized ensembles, with uniform activation and de-activation rates. The WT behaviours look similar, with some differences in propensities of phenotypes (Fig. 4(c,e)). The double-mutant on the ensembles shows a much less contrasted picture than on the Cohen's model. While Migration becomes the most likely outcome, several other phenotypes are observed, suggesting a potential variability of the effect of the double-mutation. Interestingly, even the single mutant constraints of the second ensemble are not sufficiently restrictive to guarantee the behaviour observed in Cohen's model.

Variability of propensities of phenotypes

To study the ensemble composition, we want to analyze the steady-state probabilities for each model. Depending on the results we might have a lot of visited states, which bring a dimensionality issue. We choose to represent these results using Principal Component Analysis (PCA) [START_REF] Wold | Principal component analysis[END_REF], which allows us to visualize the distribution of attractor's proportions in a reduced number of dimensions.

We apply PCA to the probability distribution of each model within the ensemble from WT and single mutants constraints, allowing us to represent their respective probability distributions (Fig. 5). The first component, representing 56% of the observed variance, shows a negative correlation between apoptotic and EMT phenotypes. The second component, representing 24% of observed variance, shows a negative correlation between EMT without Migration and EMT with Migration. The distribution of the ensemble's probability distributions is diverse, illustrating the performance of the all possible model diversity sampling. The p53 LoF/Notch GoF double mutant shows a shift towards EMT and/or Migration phenotypes, away from Apoptotic phenotypes. The alignment of models on the top-left corresponds to models which don't show any apoptotic phenotypes (∼ 96% of the models).

Conclusion

The synthesis of BNs from network architecture and dynamical constraints can lead to a multitude of admissible solutions.

In this work, we employed Answer-Set Programming to sample ensembles of diverse BNs, all possessing the same network architecture and satisfying the same set of dynamical constraints. We significantly extend the previously described methodology with the new type of biologically relevant universal constraints.

Our synthesis framework enables specifying existence and absence of reachability properties between (partial) configurations of the BN, existence of fixed points and cyclic attractors matching with observations, and universal properties on the fixed points and reachable fixed points; all these properties can be parametrized by mutation settings.

The dynamics of ensembles is explored by stochastic simulations using the new Ensemble MaBoSS simulator, which is introduced here for the first time. The ensemble-based simulations are used for computing and comparing propensities of reachable attractors under different mutations or their combinations. The result of an ensemble-based simulation represents a multidimensional distribution of the vectors of attractor probabilities, which can be aggregated by computing its mean point. Moreover, the multi-variate variance of the distribution can be explored, e.g. by applying the standard machine learning methods such as Principal Component Analysis, which can lead to the insights about the diversity of possible modelling scenarios compatible with available biological knowledge and the experimental data.

As illustrated on a biological case study, BN ensemble modelling brings an insight into the potential variability of predictions subject to model uncertainty.

In future work, we plan to address the evaluation of the diversity of sampled ensembles, with metrics helping estimate required sample size and compare sampling heuristics.

1 3 2 Fig. 1 .

 21 Fig. 1. Example of Boolean network f and its influence graph G(f ) where positive edges are with normal tip and negative edges are with bar tip.

Definition 2 (

 2 Fully-asynchronous semantics).

  We write ρ f a1 (x) := {y ∈ B n | x f -→ a1 * y} the set of configurations in transitive relation with x, with f -→ a1 * the reflexive and transitive closure of f -→ a1 .

Fig. 2 .

 2 Fig. 2. Examples of smallest K-closed hypercubes containing the configuration 100 for the BN f of dimension 3 defined by f1(x) := 1, f2(x) := x1, f3(x) := x1 ∧ ¬x3. Configurations belonging to the hypercube are in bold; these verifying the MP reachability property are boxed. The hypercube 11 * is the only one which is closed by f and minimal.

-

  an influence graph G = {[n], E + , E -) p partial observations o 1 , . . . , o p sets FP and UFP of indices of observations sets PR and URFP of couples of indices of observations: URFP ⊆ [p] 2

Fig. 3 .

 3 Fig. 3. Influence graph of Cohen's model relating 32 nodes with 159 edges, where positive edges are in green and negative in red.

Fig. 4 .

 4 Fig. 4. Simulations results for phenotypes propensities in Cohen's model (a,b), ensemble from WT constraints (c,d), and ensemble from WT and single mutants constraints (e,f), in wild-type condition (a,c,e) and double-mutant p53 LoF/Notch GoF (b,d,f).

Fig. 5 .

 5 Fig. 5. PCA representation of the steady-state distribution of each model of the ensemble from WT and single mutants constraints. Each point represents the result of one model simulation (blue one are from WT simulations, orange one from p53 LoF/Notch GoF). Large blue and orange circles highlight the position of the original single Cohen's model simulation. The triangular pattern of the distribution comes from the fact that the phenotype probabilities are located in the n-dimensional simplex.

https://maboss.curie.fr, https://github.com/colomoto/pyMaBoSS

Code, data, and notebooks at https://doi.org/10.5281/zenodo.3938904; Synthesis has been performed on 36-cores CPUs @ 2.6Ghz with 192Go of RAM; first ensemble was generated at a rate of 5s/model/CPU; second ensemble was generated at a rate of 3min/model/CPU

Acknowledgements This work has been partially supported by Agence Nationale de la Recherche in the program Investissements d'Avenir (project No. ANR-19-P3IA-0001; PRAIRIE 3IA Institute), by ANR-FNR project "AlgoRe-Cell" (ANR-16-CE12-0034), by ITMO Cancer, and by the Ministry of Science and Higher Education of the Russian Federation (project No. 14.Y26.31.0022). Experiments were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr).


Thomas Martzel

Julien Annibaletto

Pierre Millet

Etienne Pair

Morgane Sanselme

Sylvain Oudeyer

Vincent Levacher

Jean-François Brière

Organocatalytic Multicomponent Synthesis of α/β-Dipeptide Derivatives

Keywords: amino acid, asymmetric synthesis, isoxazolidinone, Meldrum's acid, organocatalysis

A straightforward multicomponent Knoevenagel-aza-Michael-Cyclocondensation reaction involving readily available hydroxamic acid-derived from naturally occurring α-amino acids allows a diversity-oriented synthesis of novel isoxazolidin-5-ones possessing an N-protected α-amino acid pendant with good to high diastereoselectivities thanks to a match effect with a chiral organocatalyst. These diversely substituted heterocycles, easily isolated as a single diastereoisomer, proved to be versatile platforms for the formation of an array of α/β-dipeptide fragments.

 The synthesis of the α/β-dipeptide motif essentially relies on the construction of the amide bond by means of coupling reagents, essentially stoichiometric, between suitably (orthogonal) protected and pre-elaborated α-amino acids (αAA) and β-amino

acids (βAA) derivatives (Scheme 1A). [START_REF] De Figueiredo | For emerging coupling reagents for amide bond formation[END_REF] Bode and co-workers have made a major achievement in this field, by developing the αketoacid-hydroxylamine (KAHA) ligation between two presynthesized AA derivatives, one bearing a α-keto acid pendant and another possessing an hydroxylamine N-terminus functional group. [7] Recently, Takemoto expended the portfolio of hydroxylamine-derived βAA to aspartic acid derivatives thanks to an organocatalytic aza-Michael addition to fumaric monoacid electrophiles. [8] This elegant strategy, making use of hydroxylamine functions as key building blocks for the elaboration of βAA, was applied to the elaboration of α/β-dipeptides albeit with moderate diastereoselective ratio (dr).

Scheme 1. A new entry to β-amino acid derivatives.

Besides these convergent sequences, we tackled an alternative divergent approach in order to populate the chemical space within the biorelevant α/β-dipeptide series (Scheme 1B). We reasoned that the multicomponent Knoevenagel-aza-Michael-Cyclocondensation (KMC) reaction would take place between readily available hydroxamic acids 1, [START_REF] Zhu | Multicomponent Reactions in Organic Synthesis[END_REF] derived from naturally occurring and enantiopure αAA, various aldehydes 2 and Meldrum's acid (MA) as a C2-synthon (first point of diversity). [START_REF] Berini | For early discovery of the MCR to form racemic isoxazolidin-5-ones from achiral aldehydes and N-oxycarbonyl-hydroxylamines[END_REF], [START_REF] Dumas | For reviews on catalytic transformations of alkylidene Meldrum's acids[END_REF] Then, this strategy would afford a straightforward elaboration of novel isoxazolidin-5-ones 3, [START_REF]For a review on catalytic enantioselective synthese of isoxazolidin-5-ones[END_REF] as masked βAA with an αAA side chain, providing eventually versatile platforms to elaborate a library of α/β-dipeptides (second point of diversity). [START_REF] Postikova | For rare examples of catalytic enantioselective synthesis (dicomponant) of -substituted N-EWG isoxazolidin-5-ones[END_REF], [START_REF] Ibrahem | For dicomponant but two steps catalytic synthesis of -substituted N-EWG isoxazolidin-5-ones[END_REF] In order to be synthetically useful, this sequence should demonstrate a significant functional group compatibility and furnish one stereoisomer. Thanks to the marked electrophilic E-mail: jean-francois.briere@insa-rouen.fr [b] Dr. M. Sanselme Laboratoire SMS -EA3233, Normandie Univ-University of Rouen, France. [c] These authors contributed equally to this work.

Supporting information for this article is given via a link at the end of the document.((Please delete this text if not appropriate)) reactivity of alkylidene MA intermediates 4, [START_REF] Dumas | For reviews on catalytic transformations of alkylidene Meldrum's acids[END_REF], [START_REF] Kaumanns | Benzylidene MA is 10 11 time more electrophilic than benzylidene malonate[END_REF] we postulated that (1) a facile domino addition process of hydroxamic acids 1 would occur and (2) the C-N bond formation between 1 and 4 would be under the influence of both a chiral nucleophile 1 and a Brønsted base organocatalysts to favor the formation one stereoisomer (Scheme 1B). We are delighted to report on the unprecedented diastereoselective multicomponent reaction (MCR) giving rise to versatile isoxazolidin-5-ones 3 thanks to the match influence of a suited organocatalyst, en route to the elaboration of various original α/β-dipeptide derivatives.

Table 1. Proof of principle and optimization. [a] Entry Deviation from the standard conditions Yield [%] dr [b] 1 Without catalyst phthalimidyl (5a) 20 (55) [c] 57 [a] Reaction conditions: 1a (0.1 mmol), PhCHO 2a (1 equiv), MA (1.5 equiv) in PhMe/CH2Cl2 (3:1, 0.025 M) at 30 °C for 5 h. Yields of both diastereoisomers determined by 1 H NMR with an internal standard. [b] dr determined by 1 H NMR on the crude product. [c] Bimolecular reaction from benzylidene Meldrum's acid 4a. [d] MA (1 equiv), (DHQ)2PHAL (20 mol%).

At the onset, it was shown that the model MCR involving the N-Boc alanine hydroxamic acid 1a (Boc-αAlaNHOH), MA and benzaldehyde 2a furnished after 5 hours the corresponding isoxazolidinone 3a in 20% yield (determined by 1 H NMR) as 57:43 mixture of diastereoisomers namely with virtually no stereoinduction (entry 1, Table 1). By performing the dicomponent reaction with benzylidene Meldrum's acid 4a, the putative intermediate involved into the formal domino aza-Michael-cyclocondensation process (see Scheme 1), similar poor dr was obtained with a better 55% yield (entry 1). The facility with which this MCR takes place is worthy of note, likely thanks to the high electrophilic reactivity of benzylidene MA, [START_REF] Kaumanns | Benzylidene MA is 10 11 time more electrophilic than benzylidene malonate[END_REF] but the ability of a catalyst to accelerate this process and promote a stereoselective sequence remains an issue. [START_REF] Berini | For early discovery of the MCR to form racemic isoxazolidin-5-ones from achiral aldehydes and N-oxycarbonyl-hydroxylamines[END_REF] To our delight the use of 10 mol% of a Brønsted base like the achiral quinuclidine allowed to improve significantly the yield to 86% and dr (71:29, entry 2). A large screening of organocatalysts (see SI) showed that the commercially available and chiral amine (DHQ)2PHAL turned out to be the more efficient catalyst leading to good 86:14 dr in 75% yield in 5 hours (entry 3 versus entries 4-7). Early optimization endeavors (see SI) showed that less polar solvents like CH2Cl2 (entries 8-9) benefited to the level of stereoinduction, while a mixture of PhMe/CH2Cl2 was required for solubility issue (entries 10-11). Noteworthy, isoxazolidone 3a proved to be somewhat unstable during the purification by silica-gel column chromatography but also in the presence of the Brønsted base in solution (vide infra). Then, some decomposition events occur during longer or forcing reaction conditions leading to eratic outcomes. However, the use of more diluted conditions and a slight excess of MA prevent the decomposition of product 3a (thanks to its high acidity, pKa = 4.8 in H2O) [START_REF] Dumas | For reviews on catalytic transformations of alkylidene Meldrum's acids[END_REF], [START_REF] Sohail | For a recent example describing the use of a 1,3-diketone that buffer nonaqueous Solutions[END_REF] even at 30 °C, a temperature allowing both a faster process and a decrease of the amount of catalyst from 20 to 10 mol% while giving improved yield of 75% and the same dr in 5 hours (entries 3 vs 11, Table 1). This shows the key role of MA to secure soft reaction conditions. Interestingly, it was observed a match/mismatch effect between chiral Cinchona-derived organocatalysts and (S)-Boc-αAlaNHOH 1a (Table 1, entries 3-4 and 6-7) which demonstrated that quinine versus quinidine derivatives were more competent to give (S,S)isoxazolidinone 3a with a good dr. Accordingly, it was proven that quinidine derived (DHQD)2PHAL catalyst was able to transform the (R)-Boc-αAlaNHOH enantiomer 1a' into isoxazolidinone (R,R)-3a' in 71% yield and 82:18 dr contrary to (DHQ)2PHAL catalyst giving 3a' only in 58% yield and 77:23 dr (Scheme 2). Eventually, it was validated that this MCR took place with no racemization and provided the major isoxazolidones 3a and 3a' as a single enantiomer (>99:1 er, see SI).

Scheme 2. Molecular diversity towards a versatile synthesis of 3a' -the other enantiomer of isoxazolidinone 3a. Yields of both diastereoisomers determined by 1 H NMR with an internal standard. [a] 47% isolated yield after column chromatography.

The scope and limitation of this novel MC-strategy was next addressed (Table 2). Pleasingly, carrying out the reaction for 24 hours on a larger scale secured the complete transformation of Boc-αAlaNHOH 1a into isoxazolidinone 3a in 88:12 dr (99% of NMR yield for both diastereoisomers). Importantly, the two diastereoisomers of isoxazolidinone 3a could be separated by flash column chromatography, as a general behavior in this series, leading to the major (S,S)-3a in 71% isolated yield (slight drop of the theoretical yield of 88%). [START_REF]diastereoisomers of isoxazolidinones 3 is a general rule in this series[END_REF] Starting from isovaleraldehyde 2b, a facile synthesis of the corresponding isoxazolidinone 3b, having βhLeu backbone (homologated β 3 -Leucine), occurred in excellent dr > 95:5 and 96% yield for the pure major stereoisomer. As a rule of thumb, in comparison to the use of aromatic aldehydes such as 2a, the isoxazolidinones 3 derived from aliphatic aldehydes 2 turned out to be more stable which allows (1) an easy purification by chromatography and (2) to carry out the reaction in higher concentration. [START_REF]In case of aliphatic aldehydes such valeraldehyde 2b, the catalyzed reaction (optimized conditions) could be performed with more conditions with only[END_REF] The MCR was applied successfully to N-Cbz-αAlaNHOH 1b (3c-3d, 64-85%, 83:17-89:11 dr) [19] and N-Fmoc-αAlaNHOH 1c (3e, 86%, 93:7 dr) to provide the corresponding 3c-3e with good to excellent stereoselectivities and yields. Then, various βhLeu-derived isoxazolidinones were constructed from isovaleraldehyde and Boc-protected proteinogenic αAA such as αPhe (3f, 71%, 89:11 dr), αVal (3f, 93%, >95:5 dr) and αPro (3i, 63%, 75:25 dr) with yields ranging from 63% to 93%. It was also shown that N-Boc glycine hydroxamic acid afforded the corresponding isoxazolidinone 3h in excellent 93% yield but low 58:42 er. In line with previous observations, this highlights a moderate capability of (DHQ)2PHAL catalyst to promote an enantioselective process from an achiral αAANHOH. [10a], [START_REF] Pair | For a rare example of enantioselective catalytic C-N bond formation to alkylidene Meldrum's acids[END_REF] Table 2. Scope and limitations. [a] [a] Reaction conditions: 1 (0.3 mmol), RCHO 2 (1 equiv), MA (1.5 equiv) in PhMe/CH2Cl2 (3:1, 0.025 M) at 30 °C for 24 h. Isolated yields of the pure major diastereoisomer after column chromatography. In parentheses, yields of both diastereoisomers determined by 1 H NMR with an internal standard. dr determined by 1 H NMR on the crude product. [b] 5 mol% of catalyst on 1 mmol scale gives 3b (89%, 94:6 dr), 3d (86%, 89:11), 3e (80%, 94:6 dr), 3r (47%, 94:6 dr). [c] 20 mol% of catalyst. [d] Absolute configuration was not determined for 3h.

An array of both aromatic aldehydes (3j-3l, 51-68%, 88:12 to 92:8 dr), linear and cyclic aliphatic aldehydes (3m-3p, 69-91%, 85:15 to > 95:5 dr) proved to be compatible with this multicomponent KaMC process, encompassing the construction of βhPhe 3n and βhLeu 3o precursors. Isoxazolidinones with ether 3q (βhSer, 61%, 83:17 dr), thioether 3r (βhMet, 41%, >95:5 dr) and NHBoc 3s (βhLys, 44%, 89:11 dr) pendants could be isolated in diastereomeric pure form albeit in slightly lower yields even with 20 mol% of catalyst, due to the instability of these aldehydes that results in an incomplete reaction. Although the homologous tryptophane derivative 3t (βhTrp) could be constructed with excellent > 95:5 dr, a low conversion was observed leading to a moderate 19% isolated yield. Worthy of note, several MCR could be carried out on 1 mmol scale to give the isoxazolidinones 3b, 3d, 3e and 3r with rather similar yields and dr by means of only 5 mol% of (DHQ)2PHAL as organocatalyst. Having an array of novel isoxazolidinones 3 in hand, decorated with various functionalities, the transformations of 3 into α/β-dipeptide derivatives was undertaken in order to probe the versatility of these masked βAA platforms (Scheme 3). By means of adapted palladium catalyzed N-O bond hydrogenolysis, followed by esterification of the obtained acids 5, in order to facilitate the purification step, several syn-α/β-dipeptides-OMe 6a-6c (86-93%), [19] together with NHBoc 6f (Boc-αAla-Boc-βhLys, 80%) and SMe 6g (Boc-αAla-Boc-βhMet, 67%) derivatives were easily synthesized. This straightforward strategy was also applied to the formation of (R,R)-dipeptide 6a' in 92% yield. [19] On the other hand, by means of orthogonal Zn/AcOH deprotection conditions, the elaboration of N-Cbz and N-Fmoc αAla-βhLeu dipeptides 6d (48%) and 6e (74%) respectively was allowed. [10b] Product 6h, displaying a sensitive but useful terminal alkyne moiety, was accessible in 50% yield under these reductive conditions (Zn/AcOH) and successfully transformed into triazole 10 in 82% yield, thanks to the copper-catalyzed 1,3-dipolar cycloaddition protocol. The readily available Boc-αAla-βhLeu carboxylic acid 5b (used as a crude product) was subjected to peptide-coupling reagents to give rise to the formation of either thioester 7 (94%), suited for further native ligation processes, or α/β/α-tripeptide 8 in 91% yield. Interestingly, taking advantage of the reactivity of the N-EWG isoxazolidinone framework, [12],[7b, 7c] the direct ring opening event of 3b with tert-butyl-alanine furnished the α/β/α-tripeptide 9, similar to 8 but with a N-OH functionality known to be a valuable moiety to complex metal in bio-environment. [21] In summary, we have developed a multicomponent synthesis of novel isoxazolidin-5-ones possessing an N-α-amino acid. Thanks to the high reactivity of transient alkylidene Meldrum's acids, a smooth and diastereoselective domino addition reaction takes place upon the match influence of a commercially available quinine derived organocatalyst allowing to achieve good to excellent dr. The corresponding isoxazolidin-5-ones, easily obtained as single diastereoisomer after purification, proved to be versatile platforms for the diversity-oriented synthesis of an array of α/β-dipeptides, as useful fragment in medicinal chemistry.

  . Martzel, Dr. J. Anniballeto, P. Millet, Dr. E. Pair, Dr. S. Oudeyer, Dr. V. Levacher and Dr. J.-F. Brière Normandie Univ, UNIROUEN, INSA Rouen, CNRS, COBRA, 76000 Rouen, France.

Scheme 3 .

 3 Scheme 3. Chemical transformations into α/β-dipeptides. Reaction conditions: [a] H2, Pd/C (1 atm), i-PrOH, 30 °C for 6a-b;15 atm, RT for 6c, 1 atm, 50 °C for 6f and TMSCHN2, MeOH, RT. [b] H2, Pd(OH)2/C (20 atm), i-PrOH/EtOAc (2:1), 60 °C for 6g. [c] Zn (40-80 equiv), THF/H2O/AcOH, 40 °C, and TMSCHN2, MeOH, RT. [d] From 3b, HCl•AlaOt-Bu, iPr2NEt, DMAP (0.5 equiv), DMF, RT. [e] From 5b, HCl•AlaOt-Bu, i-Pr2NEt, EDC•HCl, HOBt, CH2Cl2, RT. [f] PhSH, DCC, HOBt, CH2Cl2, RT. [g] BnN3, Na•Ascorbate, CuSO4•5H2O (20 mol%), t-BuOH/H2O, RT.

Acknowledgements

This work has been partially supported by INSA Rouen Normandy, University of Rouen Normandy, the Centre National de la Recherche Scientifique (CNRS), EFRD, and Labex SynOrg (ANR-11-LABX-0029), and by Region Normandie (CRUNCh network). This research was also supported by the French National Research Agency (ANR) as part of the ANR-16-CE07-0011-01 project OMaChem.

COMMUNICATION

Entry for the Table of Contents (Please choose one layout)

COMMUNICATION

The diversity: A straightforward organocatalyzed multicomponent reaction (MCR) involving readily available hydroxamic acids-derived from naturally occurring α-amino acids (αAA) allows a stereoselective and diversity-oriented synthesis of novel isoxazolidin-5-ones as versatile platforms for the elaboration of α/β-dipeptide fragments.

Thomas


Ronan Le Goff

Olivier Mahé

Ronan Le Coz-Botrel

Sylvie Malo

Jean-Michel Goupil

Jean-François Brière

Isabelle Dez

Ronan Le Goff

Insight in chitosan aerogels derivatives -Application in catalysis

Keywords: Chitosan aerogel, Reductive amination, Textural properties, Heterogeneous catalysis 1-Introduction

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

functional groups were covalently introduced and an original application the palladiumcatalyzed allylic substitution was allowed thanks to a phosphine-substituted material.

2-Materials and methods

Materials

Et2O, DCM, THF and MeCN were purified by a Pure Solv-MD 5 Solvent Purification System (activated alumina column containing a copper catalyst and molecular sieves). Anhydrous THF and MeCN were degassed (freeze-pump) before use. Absolute EtOH was used without further purification for Schiff Base formation and for chitosan and chitosan derivatives aerogel beads formation. MeOH was dried over 3Å molecular sieves when implied as solvent in the reduction step 37. 1,2-DCE was purchased in 99% grade from Alfa-Aesar and dried over 3Å molecular sieves before use. Chitosan middle-viscous (Mv = 336,000 g.mol -1 ; DA = 20%) was purchased from Sigma-Aldrich. Commercially available reagents p-tolualdehyde (Janssen; 97%), methyl 4-formylbenzoate (Alfa-Aesar; 99%), p-anisaldehyde (Fisher-Acros; 99%), 4-(4,4,5,5tetramethyl-1,3,2-dioxaborolan-2-yl)-benzaldehyde (Fisher-Acros; 97%), 1,3-dibromobenzene (Fisher-Acros; 97%), butyl lithium (Sigma-Aldrich; 1.6M solution in hexane), sodium triacetoxyborohydride (Fisher-Acros; 97%), sodium cyanoborohydride (Fisher-Acros; 95%), tetrabutylammonium borohydride (Sigma-Aldrich; 98%), 4-(dimethylamino)pyridine (DMAP) (Sigma-Aldrich; ≥99%), bromobenzene (Sigma-Aldrich; 99%), magnesium (Sigma-Aldrich; 99%), iodine (Fisher-Acros; 99%), acetic anhydride (Sigma-Aldrich; 99%), N,Obis(trimethylsilyl)acetamide (Fisher-Acros; 95%), potassium acetate (Sigma-Aldrich; ≥99%), Allylpalladium(II) chloride dimer (Sigma-Aldrich; 98%) were used without further purification. Chlorodiphenylphosphine, dimethylformamide, cinnamaldehyde, trimethylamine, dimethyl malonate 5 were distilled before use.

Thin layer chromatography (TLC) was performed on silica gel 60-F254 plates (0.1 mm) from Merck with UV detection. Chromatographic purifications were conducting using Merck silica gel Si 60 (40-63 μm). CO2 supercritical drying process was carried out with a Critical Point Dryer Polaron E3100. NMR spectra were recorded on a Bruker Avance III 400 at room temperature in CDCl3 for organic compounds and on a Bruker Avance III 500 at 70°C in a D2O/DCl mixture for chitosan derivatives. 1 H and 13 C NMR chemical shifts (δ) are reported in ppm using the residual peak of chloroform-d (7.26 for 1 H NMR; 77.16 ppm for 13 C NMR) or water (4.79 for 1 H NMR). Coupling constants (J) are reported in Hertz (Hz). The following abbreviations were used to explain the multiplicities: s = singlet, d = doublet, t = triplet, q = quartet, m = multiplet, bs = broad singlet. Solid-state MAS NMR (Magic Angle Spinning Nuclear Magnetic Resonance) spectroscopy has been performed on a Bruker Avance III spectrometer using a 4 mm double resonance probe. The spectra has been obtained using 90° single-pulse with a duration of 2.26 μs, a recycling delay of 5 s and a spinning rate of 9 kHz.

InfraRed spectra were recorded on a Perkin Elmer Spectrum One FT-IR spectrometer equipped with an ATR device, and only the strongest or structurally most important peaks are listed in cm -1 . High Resolution Mass Spectrometry was performed on a Q-TOF Xevo G2-XS WATERS spectrometer. Surfaces areas were determined by N2 physisorption at 77K using an ASAP 2020 Micromeritics apparatus. The powdered sample (200 mg) was evacuated overnight at 353K prior N2 adsorption. The total porous volume Vpore is dertermined from the volume adsorbed under relative pressure of 0.99. High resolution X-ray powder diffraction data were collected using a D8 Advance Vario 1 Brucker diffractometer (θ-2θ Bragg-Brentano mode) equipped with a Cu X-ray tube (λ = 1.540598 Å), a Ge ( 111 The absolute determination of weight and number average molar masses (respectively, Mw and Mn) of chitosan and Schiff base chitosan aerogels 1 and 2a were performed by coupling online a size exclusion chromatography (SEC), a multi-angle laser light scattering (MALLS) and a differential refractive index detector (DRI). Acetic acid/acetate buffer (pH 4.7) and LiNO3 0.1 M used as carrier, was eluted at 0.5 mL.min -1 flow rate (Flom HPLC pump 301). The SEC-MALLS technic has been described elsewhere 39. The concentration of each eluted fraction have been determined with the DRI (Shimadzu RID-6A, Japan) according to the known value of dn/dC (0.19 mL.g -1 ).

Preparation of chitosan aerogel beads

1g of chitosan was suspended in 100 mL of a solution of acetic acid (0.055 M, 330 µL of AcOH in nearly 99.67 mL of ultrapure water) and stirred for 16 h at room temperature. The chitosan solution was then filtered on a glass filtering funnel and filtrate was standing for 30 minutes until no more bubbles remains. The filtrate was then dropped into a 1 liter sodium hydroxide solution (10 g of NaOH in 1 L of ultrapure water) through a 0.8 mm gauge syringe needle and the suspension was slowly stirred for another 3 hours. Chitosan beads were then dehydrated by successive immersion in a series of water-ethanol baths of increasing alcohol concentration (10, 30, 50, 70, 90 and 100%; 15 min in each bath). Chitosan beads were finally dried under supercritical CO2 conditions to give the corresponding aerogel material.

Typical procedure for the functionalisation of chitosan aerogel beads.

Chitosan aerogel beads 1 (25 mg, 0.127 mmol, 1 eq) were suspended in 1.2 mL of an ethanol solution of aldehyde (0.508 mmol for 4 eq or 0.254 mmol for 2 eq of aldehydes a-e) and heated at 60°C without stirring. After 24h, the suspension was cooled to room temperature and the solvent was removed. Chitosan Schiff base beads, recovered quantitatively, were washed 5 times by successive immersion in 1.2 mL of an ethanol baths (15 min of orbital shaking at 750 rpm for each bath). The beads were finally dried under supercritical CO2 conditions.

The beads of Chitosan Schiff base 2a-2e (5 mg, 1 eq) were placed into a 1.2 mL vial. The vial was purged with argon for 10 min before the addition of 400 µL of anhydrous methanol. A very low amount of phenol red indicator (it gives a yellow coloration for pH below 6.8 and a red one above 8.2) was added and the suspension was slowly stirred for 2 min. DMAP was added to the yellow suspension until its coloration turned to red. Finally the pH was adjusted (turned to yellow coloration) with a small amount of acetic acid. NaBH3CN (6 eq) was added to the suspension and the vial was slowly stirred at room temperature. Addition of another amount of acetic acid could be necessary to adjust pH during the reaction (coloration turning to orange/red). After 48h, solvent was removed and beads were washed 5 times by successive immersion in 1.2 mL ethanol bath (15 min of orbital shaking at 750 rpm for each bath). The chitosan amine derivatives beads 3a-3e, recovered quantitatively, were finally dried under supercritical CO2 conditions.

General Procedure for Pd-Catalysed Allylic Substitution of (E)-1,3-Diphenyl-3acetoxyprop-1-ene with Dimethyl Malonate or morpholine 40

Chitosan phosphine Schiff base beads 2e (8.7 mg, 0.0215 mmol, 0.043 eq) and [PdCl(C3H5)]2

(2.7 mg, 0.0075 mmol, 0.015 eq) were placed into a Schlenk. The Schlenk was purged three times (vacuum then argon) then 500 µL of anhydrous and degassed MeCN were added. The suspension was stirred slowly at room temperature for 1h. KOAc (0.5 mg, 0.005 mmol, 0.01 eq), BSA (370 µL, 1.5 mmol, 3 eq), a solution of (E)-1,3-diphenylallyl acetate 4 (126 mg, 0.5 mmol, 1 eq) in 700 µL of anhydrous and degassed MeCN and dimethyl malonate 5 (1.5 mmol, 3 eq) as nucleophile were successively added. Then the suspension was slowly stirred at room temperature for 24h. Solvent was removed and beads were washed three times with 900µL of anhydrous and degassed MeCN (15 min of stirring for each). Solutions were gathered and concentrated under vacuum. The crude mixture was purified by flash chromatography (SiO2; 9:1 pentane/AcOEt) to yield 6 as colorless oil (118 mg, 73%).

Charaterisation of chitosan aerogel beads functionalisation.

The determination of substitution degree (DS) for both Schiff Base formation and reduction steps have been performed by NMR analysis in a mixture D2O/DCl (78 µL of 20 w% solution of DCl + 920 µL of D2O) at 70°C. In the case of Schiff bases, DSimine were calculated using integrals of aromatic protons from free aldehydes (obtained by hydrolysis of imine moieties during NMR samples preparation) and anomeric proton from glucosamine units. In the case of secondary amines, DSamine were calculated using integrals of both anomeric protons of glucosamine and N-benzyl glucosamine units.

Chitosan Schiff base 2a

DSimine = 92-99%. IR (neat): ν = 1636 (C=N), 1065, 1023, 812 (C-HAr) cm -1 .

Chitosan secondary amine 3a

DSamine = 81-82%. IR (neat): ν = 1065, 1022, 806 (C-HAr) cm -1 . 1 H NMR (500MHz, D2O): δ = 7.80 (d, J = 7.6 Hz, HAr), 7.74 (d, J = 7.6 Hz, HAr), 5.69-5.34 (m, Hanomeric secondary amine), 5.30-5.22 (m, Hanomeric glucosamine), 4. 91-4.75 (m), 4.56-4.47 (m), 4.38-4.26 (m), 4.23-4.13 (m), 3.74-3.57 (m), 2.76 (s, Me), 2.46 (s, MeAc).

Chitosan Schiff base 2b

Chitosan secondary amine 3c

DSamine = 83-87%. IR (neat): ν = 1514, 1247 (C-O), 1064, 1027, 816 (C-HAr) cm -1 . 1 H NMR (500MHz, D2O): δ = 7.79 (d, J = 7.6 Hz, HAr), 7.41 (d, J = 7.6 Hz, HAr), 5.49-5.29 (m, Hanomeric secondary amine), 5.28-5.17 (m, Hanomeric glucosamine), 5.02-4.62 (m), 4. 49-4.36 (m), 4.33-4.20 (m), 4.17 (s, OMe), 4.15-4.00 (m), 3.68-3.50 (m), 2.39 (s, MeAc).

Chitosan

Results and discussion

Functionalisation of chitosan aerogels

The heterogeneous chitosan modification was performed from the previously prepared chitosan aerogels 1 (vide infra) as starting material, and was carried out in a two-steps reductive amination reaction (Scheme 1). In line with Quignard pioneering investigation on chitosan aerogels 35, the imine formation took place uneventfully in EtOH given that an excess of p-tolualdehyde a (4 equivalents) was used at 60°C in order to secure both a reasonable reaction rate and a high substitution degree (DSimine) of 95% determined by NMR (Figure 1). With respect to FTIR spectra, the formation of imino chitosan 2a was confirmed by the presence of the imino band at 1636 cm -1 , and the characteristic bands for the aromatic region at 1065, 1023, 812 cm -1 . Then, the heterogeneous reduction of imine derivatives 2a into benzylated-product 3a was undertaken (Table 1).

Table 1: Conditions of reduction of Schiff base 2a

Entry

Reducing agent (eq.) Additive (eq.) Solvent NMR spectra of 3a (Figure 1) demonstrate the efficacy of the NaBH3CN protocol in MeOH.

Importantly, Cabrera-Escribano and co-workers described an interesting reduction of iminederived chitosan as suspension powder in the presence NaBH3CN and AcOH, but with a lower substitution ratio, showing the advantage of using aerogel materials 39 36. Furthermore, the resulting aerogel materials, either in the form of imine or amine (2a-3a) were conveniently recovered by simple filtration. However, the impact on the texture of the functionalized material was an issue to be questioned. So the chitosan aerogels behave as material with large surface area of 260 m 2 g -1 , build in large pores (Table 2). The textural properties data of the corresponding imines 2a are reported in Table 2 and showed that this protocol has a low impact on the textural properties of the aerogels, and no shrinkage of the structure was observed. Surface area, total porous volume and PSD (centered nearby 23 nm) are highly preserved (Table 2). The high-magnification image recorded on a chitosane aerogel 1 beads surface evidenced an irregular surface based on fibrous network structure (Figure S1 in the Supporting Information). After modification, the topography of chitosan Schiff base aerogels 2a is similar.

Figure 3 showed the XRD patterns of chitosan aerogel and chitosan derivatives aerogels. Crystallinity index (CrI) reflecting the regularity of the aggregation state of polymer chains, values were listed in Table 2. The CrI of the chitosan aerogel 1 (40%) was lower than that of 2a, showing the influence of the aromatic side groups leading to a restructuration of chitosan chains to reach an aggregate architecture with high CrI.

The starting chitosan and chitosan Schiff base aerogels 1 and 2a were used for the HPLC/SEC studies. The samples were solubilized in the mobile phase, which was 0.2 M AcOH and 0.1 M NH4OAc to decrease aggregate formation 45. The determined are in the range of 100,000 g mol -1 for chitosan derivative aerogel 2a, being 50,000 g mol -1 that of the parent chitosan aerogel 1. Consequently, no depolymerisation of the parent chitosan under the reaction conditions used takes place.

The textural properties data of the chitosan secondary amine derivative 3a are reported in Table 2 and showed that the reduction conditions poorly impact the textural properties of chitosan aerogels. As the other derivatives, 3a aerogels are formed by a loose network with a surface area of 200 m 2 g -1 and a constant distribution of porous volume at 24 nm (Table 2, Figure 2) corresponding to solids with an important macroporosity. Such open accessible surface is crucial for the catalytic performance of the material affording easy penetration of reactants, diffusivity, transport… . No particular change of the topography of 3a was evidenced in comparison to those of 2a (Figure S1 in the Supporting Information).

The XRD spectra of 3a is similar of 2a with a cristallinity index of 55% showing the reduction step do not change the chitosan packing (Table 2, Figure 3).

Having this very promising protocol in hands, the modification strategy was extended to various aldehyde derivatives (Table 3) bearing aryl ring with electron-donating, withdrawing, boronic ester and phosphine functional groups making use of NaBH3CN as reductive agent. Whatever the electronic nature of aldehyde precursors, the corresponding products were obtained with substitution degree (98-99%) for the imine formation (2a-b) and a ratio of reduction ranging from 67 to 87% for the reduced benzylated-products (3a-b, entries 1-2).

Pleasingly, the boronic ester derivatives 2d-3d were nicely tolerated with this protocol (entry 3), which offers interesting opportunities on this versatile functional group. This shows the efficiency of these experimental conditions that also afford an easy recovery of the resulting material by simple filtration. As long as the phosphine chitosan aerogels 3e is concerned, the synthetic sequence started from 3-(diphenylphosphino)benzaldehyde (e) and the imine formation succeeded in 99% DSimine to afford 2e. Importantly, a former procedure was described elsewhere and afford the synthesis of chitosan-2(diphenylphosphino)imine with a low DS 46. Then the reduction step, performed with NaBH3CN as reductive agent, afforded a priori the secondary amino chitosan aerogels 3e as testified by the low remaining amount of starting material 2e. Unfortunatly, the insolubility of 3e in the NMR solvents did not allow the determination of the DSamine. However, the IR spectra confirmed the integrity of the phosphine group on the supposing material 3e by the band at 1434 cm -1 . 31 P MAS NMR analysis of 3e was performed and showed three large peaks at 30 ppm, 11 ppm and -5 ppm 47. The more intensive peak at -5 ppm is attributed to the phosphine group of 3e, nevertheless the experimental conditions did afford the partial oxidation of the phosphine indicated by the signal at 30 ppm. Finally, the signal at 11 ppm should be attributed to the phosphine groups of the polymer network involved in hydrogen bonds. These interactions likely explained the lack of solubility of 3e in various solvents.

Application to the palladium-catalyzed Tsuji-Trost reaction

Chitosan aerogels structure being well known for their catalytic properties, chitosan aerogels flanked by phosphine moiety with a covalent link were evaluated as a potential heterogenous ligand in palladium catalyzed process 48, 49. Then, the phosphine chitosan aerogels 2e in the presence of a palladium species was engaged into a catalytic allylic substitution, one of the most important reaction for carbon-carbon and carbon-heteroatom bond formation 50 (Scheme 2).

Scheme 2. Palladium catalyzed allylic substitution

Typically, the catalytic aerogels were prepared by mixing in acetonitrile, aerogels 2e (4.3 mol%) with the palladium source [PdCl(C3H5)]2 (1.5 mol%) for 1 hour under Ar. The allylation reaction between (E)-1,3-di-phenyl-3-acetoxyprop-1-ene 4 and dimethyl malonate 5 using BSA as base was selected as a model reaction. The reaction was completed after 24 h at 20°C, affording the expected allylated malonate 6 in good 76% yield (Table 4, entry 1). No enantiomeric excess was mesured by chiral HPLC. The reaction was also performed with chitosan aerogels 1 and 2a, but in these conditions the corresponding product 6 was not observed, indicating that the no active metallic-complex was formed in this case (Table 4, entries 5 and 6). Recycling tests were next performed. They showed a high catalytic activity of the chitosan aerogels for 3 cycles (100% conversion) and a decrease in activity is observed from the fourth run, leading only to a partial conversion of 73% even after 72h (Table 4, entry 4). Detection of palladium by ICP-MS in the crude product 6 showed a low leaching of the metal during the recycling of the catalytic material (< 3%) that did not explain the catalytic deactivation 51.

This analysis showed the high coordination of Pd to the biopolymer-supported Schiff-base ligand and provides evidence of chitosan acting as a good scaffold for Pd II catalyst precursors as already observed in the nice work of Smith and al. who applied successfully chitosan-2(diphenylphosphino)imine ligand in Suzuki-Miyaura and Heck cross-coupling reactions 46.

Conclusion

An efficient, mild and non-degradative procedure to prepare secondary amino chitosan derivatives, by reacting chitosan aerogels with arylaldehydes under heterogeneous conditions, is described. This convenient procedure was successfully applied to achieve reaction of several aldehydes with chitosan aerogels to reach high degree of substitution during the imineformation process. Then, the reduction of chitosan base Schiff derivatives succeeded by using NaBH3CN as reductive agent and high reduction ratio were obtained (>80%) to afford functional secondary amino chitosans aerogels. The results obtained in this work led to the conclusion that this heterogeneous procedure affords regioselective and quantitative chitosan modifications to get a controlled well-defined macromolecular structure.

Moreover, keeping in mind to preserve the original properties of the chitosan, this study focused on the impact of the experimental conditions on the physical and textural properties of aerogels.

Textural analyses showed that the successive heterogeneous chemical modifications of chitosan aerogels do not lead to a collapse of the polymer structure and aerogels keeps a surface of about 250 m 2 g -1 with a distribution of mesopores centered at 23 nm. XRD analyses showed that, after modification, the chitosan chains were restructured to form more regular aggregate architecture, with a high CrI.

Finally, the heterogeneous and efficient process to prepare macroporous aerogels afforded the synthesis of covalently link phosphine based chitosan aerogels, which were used as phosphorylated ligand in the Pd-catalysed allylic substitution. Such resulted Pd-catalytic opened materials allow by high diffusion properties and afford efficient reaction in mild conditions. The catalyst could be recycled 3 times.

  ) monochromator (Johansson type) and a Lynx Eye detector. XRD data were collected over a [4°-40°] 2θ range with an angular increment of 0.0105°, and a step time of 2 seconds (20 h/scan) at room temperature. The crystallinity index (CrI) of the chitosan samples was determined according to the method of Focher et al. 38 using the equation: CrI110 = (I110 -Iam) / I110 × 100 ; where I110 is the maximum intensity (2 = 20°) of the (110) lattice diffraction and Iam is the intensity of amorphous diffraction at 2 = 16° of chitosan.

  DSimine = 98-99%. IR (neat): ν = 1715 (C=O), 1643 (C=N), 1282 (C-O), 1109, 1069, 1045, 1016, 766 (C-HAr) cm -1 .Chitosan secondary amine 3bDSamine = 67-74%. IR (neat): ν = 1713 (C=O), 1282 (C-O), 1108, 1067, 1019, 763 (C-HAr) cm -1 . 1 H NMR (500MHz, D2O): δ = 8.41 (bs, HAr), 7.97 (d, J = 7.3 Hz, HAr), 5.51-5.36 (m, Hanomeric secondary amine), 5.30-5.22 (m, Hanomeric glucosamine), 5. 3.68 (s, OMe), 2.40 (s, MeAc).

  DSimine 99%. IR (neat): ν = 1605 (C=N), 1514, 1253 (C-O), 1168, 1065, 1024, 830 (C-HAr) cm -1 .

  Schiff base 2d DSimine 95%. IR (neat): ν = 1641 (C=N), 1360 (B-O), 1143 (C-O), 1082, 1067, 1019, 857 (C-HAr) cm -1 .Chitosan secondary amine 3dDSamine = 86-89%. IR (neat): ν = 1360 (B-O), 1143 (C-O), 1082, 1068, 1020, 857 (C-HAr) cm -1 . 1 H NMR (500MHz, D2O): δ = 8.17 (d, J = 6.8 Hz, HAr), 7.86 (d, J = 6.9 Hz, HAr), 5.53-5.30 (m, Hanomeric secondary amine), 5.30-5.16 (m, Hanomeric glucosamine), 5.04-4.60 (m), 4.54-3.84 (m), 3.69-3.46 (m), 2.38 (s, MeAc), 1.52 (s, Mepinacol).Chitosan Schiff base 2eDSimine = 97-98%. IR (neat): ν = 1644 (C=N), 1434 (P-CAr), 1147, 1067, 1041, 1026, 742 (C-HAr), 684 (C-HAr) cm -1 .Chitosan secondary amine 3e IR (neat): ν = 1434 (P-CAr),1147 , 1067 , 1041 , 1026, 742 (C-HAr), 742 (C-HAr), 684 (C-HAr) cm -1 .

  Scheme 1. Model reaction for chitosan aerogels modifications

a

  DSamine of the obtained Bn-amine 3a determined by 1 H NMR in DCl/D2O at 70°C. b in situ formation of (nBu)4BH(OAc)3 by mixing (nBu)4BH4 and AcOH in 1:3 ratio.Various acid-stable reducing agents were tested at room temperature for 48h. NaBH(OAc)3 furnished a better reduction ratio to reach a DSamine of 51% in the presence of acetic acid (entries 1 and 2). The more soluble (nBu)4BH(OAc)3, formed in situ by mixing (nBu)4BH4 and AcOH, allowed an excellent 81% DSamine (entries 3 and 4). Interestingly, the commercially available NaBH3CN led to 82% DSamine in MeOH without the need of acetic acid additive. This is an important achievement as long as an early hydrolysis of the imine moieties favored by AcOH was met in some experiments giving some reproducibility issues. This phenomenon was reduced by working in 1,2-dichloroethane (DCE) instead of MeOH for NaBH(OAc)3 for instance, but proved to be ineffective for NaBH3CN (entries 4 versus 5). The representative 1 H

Figure 1 .

 1 Figure 1. Liquid 1 H NMR spectra of (a) chitosan aerogels 1; (b) aerogel chitosan derivatives 2a corresponding to the hydrolysed imine product namely chitosan 1 and para-tolualdehyde a (aromatic and anomeric protons are indicated) and (c) 3a (anomeric protons are indicated).

Figure 2 :

 2 Figure 2: N2 adsorption isotherms of chitosan aerogel 1 (a) and chitosan derivatives aerogels

Figure 3 .

 3 Figure 3. XRD spectra of chitosan aerogels 1 and aerogel chitosan derivatives 2a and 3a

Table 2 :

 2 Physical

and textural properties of aerogels 1, 2a and 3a

Table 3 :

 3 Synthesis of functionalized Chitosan aerogels 2b-2e and 3b-3e

	Entry	Aldehydes	Chitosan Schiff	DSimine b	Chitosan Amine	DSamine b
			base derivatives a		derivatives c

a General conditions: Chitosan aerogel beads (25mg, 0.127mmol, 1eq), aldehyde (0.508mmol, 4eq) in 1.2mL of absolute EtOH at 60°C under thermal stirring. b DSimine and DSamine determined by 1 H NMR in DCl/D2O at 70°C. c General conditions: Schiff base 2 (5mg, 1eq), NaBH3CN (6eq) in 400µL of solvent at room temperature under flipping stirring.

Table 4 :

 4 Pd-catalyzed allylic alkylation of 1,3-diphenylpropenyl acetate 4 with dimethyl malonate 5 using ligands 1, 2a and 2e aerogels a

						(4.3 mol%)	2e
						(1.5 mol%) [PdCl(C 3 H 5 )] 2	O	O
	OAc		O	O	(1mol%) KOAc (3eq.) BSA	MeO	OMe
	Ph	Ph	+	MeO	OMe	MeCN, r.t., 24h	Ph	Ph
	4 (1eq.)			5 (3eq.)			6

Acknowledgements

The authors acknowledge the financial support of the INC3 M (FR3038), the french Agence Nationale de la Recherche (ANR), through the program "Investissements d'Avenir" (ANR-10-LABX-09-01), LabEx EMC 3 , and the region NORMANDIE.

Data Availability

The raw data required to reproduce these findings cannot be shared at this time due to technical and time limitations.