Charles Séjourné

Romain Couillet

Pierre Comon

A LARGE-DIMENSIONAL ANALYSIS OF SYMMETRIC SNE

Keywords: dimensionality reduction, high-dimension, non-convex, random matrix theory, machine learning

Stochastic Neighbour Embedding methods (SNE, t-SNE) aim at finding a faithful low-dimensional representation of a high-dimensional dataset. Despite their popularity, being solution to a non-convex optimization, the behavior of these tools is not well understood. This work provides first answers by leveraging a large dimensional statistics approach, where the number n and dimension p of the large-dimensional data are of the same magnitude. We derive and study the canonical equation verified by the critical points of this non-convex optimization problem. The study notably reveals that, in a simple setup, the achievable SNE solutions correspond to a subset of those critical points. In particular, when the clusters composing the dataset are balanced in size, these solutions are symmetrical and assume closed-form expressions.

As a major conclusion, the analysis rigorously proves a long-standing heuristic statement on the "proper normalization" of the symmetric SNE: out of two natural normalization choices, only the claimed proper one leads to non-trivial solutions.

INTRODUCTION

Stochastic Neighbor Embeddings (SNE) [START_REF] Hinton | Stochastic neighbor embedding[END_REF] and its variant t-SNE [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF] are dimensionality reduction methods used for visualization purposes (usually in 2 or 3 dimensions) as a pre-analysis tool for data scientists to get a grasp on their (usually numerous and large dimensional) data. The lowdimensional representations achieved by those algorithms stem from a non-convex optimization problem, traditionally solved via gradient-descent. As such, it is far from trivial to predict the resulting solution, let alone to interpret the associated visualization. Concretely, the following natural concern arises: considering 3 equidistant vectors in R p with p > 2, 1 how would SNE or t-SNE represent those points when projected in R? The equidistant property cannot hold in R and it is unclear which representation SNE will retain. The subject of the present article is precisely to cast a first light on this riddle.

Despite the underlying non-convex optimization from which SNE solutions arise, recent tractable analyses of SNE algorithms emerged in the machine learning literature. Using dynamical system techniques, [START_REF] Linderman | Clustering with t-sne, provably[END_REF] highlights a "shrinkage" property of the low-dimensional representations: under mild assumptions, within each cluster, the data point projections converge to each other as the gradient descent progresses. This is however not sufficient to ensure consistent clustering as different clusters are not guaranteed not to merge. The analysis of [START_REF] Arora | An analysis of the t-sne algorithm for data visualization[END_REF] then builds on [START_REF] Linderman | Clustering with t-sne, provably[END_REF] and shows that the cluster centroids in the low-dimensional representation space remain separated; combined with [START_REF] Linderman | Clustering with t-sne, provably[END_REF], this ensures that the representations form well-separated clusters.

While those articles establish the consistency of the SNE methods, they do not provide an accurate view on the solution landscape: how many solutions exist? With which associated cost? Is gradient descent guaranteed to fall in a global minimum? The present work addresses some of these aspects. This however demands to break the problems of non-linearity and non-convexity of the cost function to be minimized. To this end, we use a large n, p data assumption and exploit large dimensional asymptotics. For simplicity of exposition and for the sake of retrieving insightful results, the n data follow a mixture of k large dimensional Dirac masses in R p .

Our main results may be summarized as follows:

• we obtain a limiting vector equation as, n, p grow large at the same rate, satisfied by the critical points of the underlying SNE optimization problem, and analyze their stability;

• the limiting equation appears to slightly but fundamentally differ depending on the specific normalization of the symmetric SNE method initially introduced by [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF]: as a result, we reach a first rigorous argument in favor for the SNE normalization effectively used today in practice;

• the observed low-dimensional outputs are shown to be severely biased by the number of components of each class, thereby disrupting the original SNE objective;

• in a balanced class-size setting, the limiting vector equation assumes a closed-form expression, which reveals the surprising presence of a continuum of saddle points and, as a result, a very slow convergence of SNE, when compared to faster convergence in unbalanced scenarios.

The remainder of the article presents the symmetric SNE algorithm and the large-dimensional setup used for our analysis (Section 2), before introducing our main results under the form of the limiting vector equation satisfied by the symmetric SNE critical points (Section 3) which is analyzed in a simple setting and turned into a closed-form expression. For readability, in what follows, "SSNE" will refer to the symmetric SNE [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF].

PROBLEM FORMULATION

Let x1, . . . , xn ∈ R p be input high-dimensional data vectors, and X = [x1, . . . , xn] T ∈ R n×p . The goal of SNE is to associate with each xi ∈ R p a vector yi ∈ R d , for d < p (in general d ∈ {2, 3}, but we shall subsequently also consider d = 1 for simplicity), such that Y = [y1, . . . , yn] T ∈ R n×d preserves as much of the significant structure of X in R d as possible.

The SSNE algorithm

The idea of SSNE as devised by [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF] is to build two joint probability distributions P and Q ∈ R n×n for the data points X and their (yet to be found) low-dimensional representations Y , respectively. Starting from some initialization point, the representations Y are then iteratively tuned in order for the density Q to accurately approximate P .

Precisely, for some positive and non-increasing functions f : R → R and g : R → R, the density matrix Q is defined as

Qij = G ij k =l G kl , i = j 0 , i = j , Gij = g( yi -yj 2 ) while P is alternatively defined, for Fij = f ( 1 p xi -xj 2 ), as Pij = i =j      P MN;ij = F ij k =l F kl P RMNS;ij = 1 2n F ij l =i F il + F ij l =j F lj , Pii = 0
where MN stands for "matrix normalization" and RMNS for "row-matrix normalization and symmetrization", discussed in [START_REF] Hinton | Stochastic neighbor embedding[END_REF] as two normalization procedures for P . The authors in [START_REF] Hinton | Stochastic neighbor embedding[END_REF] empirically claim that the RMNS solution yields better visualizations, as it is likely more robust to outliers. One objective of the article is to theoretically support this statement.

In the original SNE algorithm [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF], both f and g are taken to be f (t) = g(t) = exp(-t/2), thereby mimicking a Gaussian distribution for both low and high dimensional data representations. The t-SNE approach rather relies on a Student-t modification of the method in which g(t) = 1/(1 + t), the justification of which is also heuristic rather than based on theoretical supports. For simplicity of exposition, in the sequel, we restrict ourselves to the case where g(t) = exp(-t/2).

The low-dimensional representation matrix Y is then chosen in such a way to minimize the Kullback-Leibler divergence between P and

Q C(X, Y ) = KL(P ||Q) = i,j Pij log Pij Qij . ( 1
)
Not being a convex problem in Y , a local minimum is then determined using a mere gradient descent procedure. A pleasant advantage of this formulation is that, for g(t) = exp(-t/2) as considered here, the gradient at running point Y is simply defined as [2]

∂C ∂yi = 4 n j=1 (Pij -Qij)(yi -yj), i = 1, . . . , n. ( 2
)
The solutions of the SNE problem, which correspond to the local minima of C(X, Y ), are therefore associated with some of the solutions of the set {∂C/∂yi = 0, i = 1, . . . , n}, i.e., the critical points of C. The objective of the article is to precisely characterize these critical points, under a statistically convenient setting. To this end, it is convenient to provide an expression of the Hessian matrix H(Y ) at any given point Y . Specifically, we find that

H(Y ) = 4 (Q -P ) ⊗ I d - 2 γp S - 4 γ 2 p T 1n1 T n T T -D ((Q -P ) ⊗ I d )1 nd - 2 γp S1 nd
where γp = k =l G kl is the normalization constant ensuring that Q is a probability distribution, 1n is the vector of ones of size n, ⊗ is the Kronecker product, D(v) is the diagonal matrix with vector v on the diagonal, and matrices S and T of size nd × nd and nd × n respectively defined block-wise as

{S i×d+l,j×d+m } 0≤l,m≤d-1 = Gij(yj -yi)(yj -yi) T {T i×d+l,j } 0≤l≤d-1 = Gij(yi -yj).

The large dimensional setting

As the purpose of SNE is to visualize large dimensional data of R p in the low dimensional space R d , the article considers the setting where p, n are arbitrarily large and of the same order of magnitude, while d is fixed, i.e., n, p → ∞ with p/n → γ > 0. We further assume the xi's drawn from

L = k =1 c δµ , c1, . . . , c k > 0, k =1 c = 1
that is, a mixture of k masses located at µ1, . . . , µ k ∈ R p . Moreover, as p → ∞, µ = O(1) for each . 2 We may write

X = Jµ T , µ = [µ1, . . . , µ k ] ∈ R p×k , J = [j1, . . . , j k ] ∈ N n×k
where [j ]i = δ {x i =µ } is the indicator vector of class . In particular, denoting n ≡ |{xi = µ }|, by the law of large numbers, n /n → c almost surely.

As shall be seen, this simple setup already allows one to retrieve interesting insights on the SSNE algorithm behavior.

MAIN RESULTS

A few notations need be defined before introducing our main results. Since all identical xi ∈ R p (say xi = µ ) are expected to have the same representation yi ∈ R d , Y must be of the form Y = J Ỹ , where Ỹ ∈ R k×d . We similarly denote Q ∈

R k×k the matrix such that Q = J QJ T -f (0)
γp In, where the diagonal centering is needed since Q has a null diagonal.

Critical points

With these notations and the assumptions of the previous section at hand, we have the following result.

Theorem 1 (Asymptotic null gradient condition). As n, p → ∞, every solution Ỹ of (2) satisfies

L Ỹ → 0, L ≡ D(1n • )AD(1n • ) -D(1n • )D(A1n • ) (3)
where

A = 1 k 1 T k -h(M ) -n(n -1) Q h(M ) = f (0) 2pnf (0) × 2 n 1 T n• M 1n • 1 k 1 T k -2nM, P = P MN M Dn • 1 k 1 T k + 1 k 1 T k Dn • M -2nM, P = P RMNS for 1n • = [n1, . . . , n k ] T , M = { µa -µ b 2 } k
a,b=1 ∈ R k×k and Dv is the diagonal matrix with vector v on the diagonal.

Sketch of Proof. Exploiting µ

= O(1) so that 1 p xixj 2 → 0 for all i, j, following the ideas of [START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF], a Taylor expansion on the entries of P is performed in the large n, p limit which, along with the symmetry of the problem, reduces { ∂C ∂y i = 0, 1 ≤ i ≤ n} from an n-to a k-dimensional equation with f only involved through f (0) and f (0).

It is interesting to observe that Theorem 1 holds identically under both MN and RMNS normalizations, but for the expression of the function h(M ) which contains the structural information about the data statistics µ1, . . . , µ k .

According to Theorem 1, for large n, p, determining the critical points of the SSNE problem boils down to solving the k-dimensional equation L Ỹ = 0 in which the non-linear function f appearing in P is now reduced to its derivatives f (0) and f (0). Unlike f though, since d is small, the nonlinear function g, intervening in Qij, remains encapsulated within Q and cannot be simplified. Equation L Ỹ = 0 thus remains difficult to analyze due to the involved non-linear relation between Y and Q.

Simplified setting

For readability and ease of interpretation, we consider from now on the case where k = 3, d = 1, and µa = mea, for ea ∈ R p the a-th canonical vector of R p and m > 0. By a common shift-invariance of the entries of Ỹ solution to L Ỹ = 0, Ỹ may be parametrized only with (say) ỹ1 -ỹ2 and ỹ1 -ỹ3.

Under this simplified setting, the MN approach can be shown to provide trivial solutions: Proposition 1. In the MN case, any Ỹ ∈ R 3 such that ỹi = ỹj for at least one combination of i = j ∈ {1, 2, 3} is a critical point of C( Ỹ).

Sketch of Proof.

The result follows from direct calculus, fundamentally exploiting the (here quite unfavorable) symmetrical nature of P MN (unlike for P RMNS , in this setting, P MN,ij is constant for all i = j and all i = j).

The MN critical points raised in Proposition 1 are highly undesirable as they "merge" distinct clusters of R p together in R d according to the cardinality of the classes. Figure 2, which reports the landscape of L Ỹ under both MN and RMNS settings in a unbalanced cluster-size scenario, effectively pinpoints those critical points in the MN setting and shows that a symmetrical pair of them are indeed local minima. This is quite unlike the RMNS method which does not lead to trivial solutions, suggesting that RMNS is always a preferred method. This is quite at odds with the seemingly "more natural" MN approach to minimizing the Kullbach-Liebler divergence in (1); as pointed out in the proof sketch of Proposition 1, RMNS here paradoxically benefits from its being "less symmetrical" than MN.

Figure 2 provides extra empirical information. In particular, it identifies a total of seven MN but only five RMNS critical points. Excluding the trivial ỹ1 = ỹ2 = ỹ3 solution, the other critical points in fact reduce to symmetrical pairs (changing the signs of all ỹi's in a solution brings the other). Among the two RMNS pairs, only one asymptotically corresponds to a minimum (green bullets in the figure), which we empirically observed to be associated with the scenario where, in the visualization of dimension d = 1, the central cluster ỹi is the one with lowest ci value (see Figure 2-(bottom)); this local minimum is here global.

A complete understanding of all critical points remains nonetheless tedious, unless the clusters have balanced sizes c1 = . . . = c k , which the next section is devoted to explore.

Balanced mixture (c1

= . . . = c k ) for k = 3, d = 1
In the setting of the previous section, let now c1 = c2 = c3. It is easy to observe that MN and RMNS now coincide; indeed, we have in this case 1n

• = (n/k)1 k so that 2 n 1 T n• M 1n • 1 k 1 T k = M Dn • 1 k 1 T k + 1 k 1 T k Dn • M.
What is less obvious though is that this balanced cluster size scenario triggers a strongly degenerate behavior of the solutions Ỹ .

Theorem 2. Under the above assumptions, for both MN and RMNS, the set of points Ỹ satisfying L Ỹ = 0 is the union of the singleton {ỹ1 = ỹ2 = ỹ3} and of the points parametrized by the ellipse As such, Theorem 2 states that there asymptotically is an infinity of critical points. Figure 3 corroborates the result and numerically suggests, through the local evaluation of the Hessian H(Y ), that the ellipse of critical points is the successive union of saddle-point and local-minimum contiguous regions.

ỹ1 -ỹ3 = -2f (0) 15pf (0) m 3 cos t - √ 3 sin t , ỹ1 -ỹ2 = -2f (0) 15pf (0) m 3 cos t + √ 3 sin t , t ∈ R.
The presence of contiguous critical points has a severe impact on convergence, as depicted by the red paths in Figures 23for which we use the same initialization points and the same gradient descent hyperparameters. While the isolated minima in the unbalanced cluster-size scenario is a strong attractor (Figure 2) and is ultimately reached at the end of the descent, in the balanced cluster-size scenario gradient descent first reaches an arbitrary point of the limiting ellipse (Figure 3) before slowly moving towards one of the symmetrical solutions where |ỹ1 -ỹ2| = |ỹ2 -ỹ3| (or any recombination of indices); these symmetrical solutions are symbolized in red diamonds in Figure 3 and correspond to the scenario depicted in Figure 1-(top)). After reaching the ellipse though, as n, p → ∞, gradient descent slows down significantly, to the point of barely progressing as soon as n, p are of the order of a few thousands: therefore, in this

CONCLUDING REMARKS

The article demonstrates that, despite the non-convex nature of the stochastic neighborhood embeddings algorithm, large dimensional statistics is able to capture its dominant behavior. This behavior is largely non-trivial and quite counterintuitive: (i) the most symmetrical scenario suffers from degeneracy and failure to converge efficiently, (ii) the stable R d visual outputs do not appropriately translate the distances between clusters in the original space R p but are impacted by cluster sizes, an undesirable feature, and (iii) the advocated RMNS approach, preferred over the nonetheless "more natural" MN flavor, only relies in a favorable asymmetrical behavior of the distribution matrix P of the affinities f ( 1 p xi-xj 2 ) of the data in ambient space.

The work yet only scratches the surface of a more complete understanding of SNE methods: many questions remain open as to the behavior under more realistic noisy data inputs, the extension to the utmost popular (and even more asymmetrical) t-SNE method, and as to "what real data visualizations really tell?".

Fig. 1 :

 1 Fig. 1: Visualization of SSNE solutions Ỹ = [ỹ1, ỹ2, ỹ3] T , for d = 1, k = 3, p = 1 024 and µa = mea (ea the a-th canonical vector of R p ) with m = 5. (Top) Any of the 6 symmetrical solutions in the balanced setting (n1 = n2 = n3 = 500). (Bottom) Most typical solution for the unbalanced setting with n1 = 300, n2 = 700, n3 = 500.

Fig. 2 :

 2 Fig. 2: Value of log L Ỹ 2 for MN versus RMNS; n1 = 300, n2 = 700, n3 = 500, m = 5, p = 1024. Local minima in green bullets. Saddle points in black crosses. Gradient descent sample path in red dotted line.

Fig. 3 :

 3 Fig. 3: Value of log L Ỹ 2 for n1 = n2 = n3 = 500, m = 5, p = 1024. Local minima in green bullets. Saddle points in black crosses. Gradient descent sample path in red dots. Symmetrical solutions achieved at finite n, p values in red diamonds.

This condition is crucial and mimics the non-trivial setting of machine learning classification of e.g.,[START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF] where x i ∼ N (µa, σ 2 Ip) with µa , σ 2 = O(1), but here in the limit where σ 2 = 0.

Couillet's work is partially supported by the ANR-MIAI Large-DATA chair at University Grenoble-Alpes, and the HUAWEI LarDist project.


Anna Lagunas

Christine Belloir

Loïc Briand

Pau Gorostiza

Josep Samitier

Odorant binding changes the electrical properties of olfactory receptors at the nanoscale

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

Odorant Binding Changes the Electrical Properties of Olfactory Receptors at the Nanoscale

Supporting Information Placeholder

ABSTRACT:

The transduction of odorant binding into cellular signaling by olfactory receptors (ORs) is not understood and knowing its mechanism would enable developing new pharmacology and biohybrid electronic detectors of volatile organic compounds bearing high sensitivity and selectivity. The electrical characterization of ORs in bulk experiments is subject to microscopic models and assumptions. We have directly determined the nanoscale electrical properties of ORs immobilized in a fixed orientation, and their change upon odorant binding, using electrochemical scanning tunneling microscopy (EC-STM) in near-physiological conditions. Recordings of current versus time, distance, and electrochemical potential allows determining the OR impedance parameters and their dependence with odorant binding. Our results allow validating OR structural-electrostatic models and their functional activation processes.

Olfactory receptors (ORs) comprise the largest multigene family in the vertebrates, with about 400 genes identified in humans. [START_REF] Buck | A Novel Multigene Family May Encode Odorant Receptors: A Molecular Basis for Odor Recognition[END_REF][START_REF] Malnic | The Human Olfactory Receptor Gene Family[END_REF] They are expressed primarily by olfactory sensory neurons located in the olfactory epithelium in the nasal cavity and are responsible for odorant detection. Moreover, the expression of ORs in other tissues have been reported, where they exert distinct biological functions. [START_REF] Feldmesser | Widespread Ectopic Expression of Olfactory Receptor Genes[END_REF][START_REF] Lee | Therapeutic Potential of Ectopic Olfactory and Taste Receptors[END_REF] ORs belong to the class A (rhodopsin-like) family of G protein-coupled receptors (GPCRs). Around half of the GPCRs are ORs. GPCRs are the most abundant membrane proteins having widespread, significant roles in signal transduction in cells. Therefore, they are a major pharmacological target, with approximately the 40% of approved drugs on the market targeting GPCRs. [START_REF] Hutchings | Opportunities for therapeutic antibodies directed at G-protein-coupled receptors[END_REF][START_REF] Hauser | Trends in GPCR Drug Discovery: New Agents, Targets and Indications Nat[END_REF] The mechanism of olfactory transduction in the main olfactory epithelium involves OR switching from a conformationally inactive state towards an active state upon ligand binding, which couples with the intracellular G protein Golf and activates adenylyl cyclase type-III. This leads to an increase of the intraciliary level of cyclic adenosine monophosphate (cAMP) opening cyclic-nucleotide-gated (CNG) non-selective cation channels. The subsequent ionic flux causes membrane depolarization, the activation of Na + -channels and consequently, action potential firing. [START_REF] Boccaccio | Patch-Clamp Recordings from Mouse Olfactory Sensory Neurons[END_REF] ORs show high selectivity and sensitivity towards odorant detection, 8 a characteristic that prompted the development of biohybrid sensors using ORs for the detection of volatile organic compounds. [START_REF] Yang | Nanodisc-Based Bioelectronic Nose Using Olfactory Receptor Produced in Escherichia coli for the Assessment of the Death-Associated Odor Cadaverine[END_REF][START_REF] Yamada | Highly Sensitive VOC Detectors Using Insect Olfactory Receptors Reconstituted into Lipid Bilayers[END_REF] This made the electrical characterization of ORs a subject of interest, which has been mainly addressed by electrochemical impedance spectroscopy (EIS) in bulk experiments. [START_REF] Hou | A novel detection strategy for odorant molecules based on controlled bioengineering of rat olfactory receptor I7[END_REF][START_REF] Khadka | An UltraSensitive Electrochemical Impedance-Based Biosensor Using Insect Odorant Receptors to Detect Odorants[END_REF] However, their microscopic interpretation relies on theoretical models of the protein electric properties and is subject to assumptions. These models are generally based on equivalent impedance networks trying to link the dipoles and charges of the backbone and side chains of the receptor in its native and activated state, with the experimental results obtained. [START_REF] Alfinito | The Emerging Science of Protein-based Electronic Devices[END_REF][START_REF] Alfinito | Current-Voltage Characteristics of Seven-Helix Proteins from a Cubic Array of Amino Acids[END_REF] Here, we have determined directly the nanoscale electrical properties of human OR hOR1A1 in the presence of its cognate ligand dihydrojasmone, using electrochemical scanning tunneling microscopy (EC-STM) under bipotentiostatic control. In particular, we have measured current-voltage (I-V), current-time (I-t) and current-distance (I-z) characteristics of the receptor immobilized in a fixed orientation in a near-physiological environment. ORs behave as parallel resistor-capacitor (RC) circuits with R and C decreasing by 20 % upon ligand binding. Our results will allow validating OR structural-electrostatic models and their functional activation including odorant binding at the extracellular side and guanosine triphosphate (GTP) binding and G protein release at the intracellular side. hOR1A1 was overexpressed in a stable tetracycline inducible HEK293S GnTI -cell line as previously described. [START_REF] Belloir | Biophysical and Functional Characterization of the Human Olfactory Receptor OR1A1 Expressed in a Mammalian Inducible Cell Line[END_REF] The receptor was engineered by inserting the epitope tags rho1D4 at the C-terminal and FLAG at the N-terminal to allow its purification and detection. Circular dichroism analysis demonstrated that detergent-solubilized FLAG-rho1D4-tagged hOR1A1was properly folded into α-helical structure, as expected for the secondary structures of a GPCR. [START_REF] Erdogmus | Helix 8 is the Essential Structural Motif of Mechanosensitive GPCRs[END_REF] Receptor functionality was assessed by ligand binding using an intrinsic tryptophan fluorescence assay, which revealed the receptor bound its cognate odorant, dihydrojasmone, with an affinity in the micromolar range exhibiting a Kd value of 1.4 ± 0.5 µM (Figure 1A).

We imaged by atomic force microscopy (AFM) and measured the single-protein current-bias voltage (I-V) characteristics of hOR1A1 in a physiological environment (50 mM phosphate buffer pH 7.4) by EC-STM. The experiments were conducted under bipotentiostatic control of the probe and sample electrodes versus an Ag/AgCl (SSC) reference electrode, with and without the presence of dihydrojasmone. I-V measurements were performed at a fixed separation between the probe and sample (given by a setpoint current of 0.4 nA) to avoid physical contact between the STM probe and the protein. [START_REF] Artés | Current-Voltage Characteristics and Transition Voltage Spectroscopy of Individual Redox Proteins[END_REF] The STM feedback loop allowed fixing the setpoint current and it was turned off during 0.86 s to perform each I-V recording at a rate of 581 mV/s. This avoids the biasing that might be introduced by the contact geometry in contact mode, [START_REF] Zhang | Role of Contacts in Long-Range Protein Conductance[END_REF] and provides a more flexible configuration. A uniform orientation of the receptor is essential in single molecule experiments, [START_REF] Lagunas | Long Distance Electron Transfer through the Aqueous Solution between Redox Partner Proteins[END_REF][START_REF] Ortiz | Distance and Potential Dependence of Charge Transport Through the P700 Reaction Center of Photosystem I[END_REF] and to achieve it hOR1A1 was immobilized on the Au(111) electrode of the EC-STM by half anti-Rhodopsin antibody, against the C-terminal tag Rho1D4, [START_REF] Sharma | Half Antibody Fragments Improve Biosensor Sensitivity without Loss of Selectivity[END_REF] (Figure 1B). Coverage of the substrate was monitored by atomic force microscopy imaging (Figure S1). Height histograms show a peak centered at 4.8 ± 0.9 nm after the incubation with the hOR1A1, which was attributed to the receptor, in agreement with the sequence homology data obtained for hOR1A1. [START_REF] Kelley | The Phyre2 Web Portal for Protein Modeling, Prediction and Analysis[END_REF] Sample potential (US) and initial probe potential (initial UP) were set at 0.25 V and 0.45 V, respectively. Faradaic leakage current was maintained below a few pA through probe insulation. [START_REF] Artés | Transistor-like Behavior of Single Metalloprotein Junctions[END_REF] To obtain each I-V plot, the probe was positioned over the sample at a current set point of 0.4 nA, the feedback loop disconnected and the probe current recorded while a voltage ramp is applied to the EC-STM probe. Up to 150 I-V curves were recorded at different sample positions by sweeping the probe potential back and forth from 0.15 to 0.40 V (bias (Ubias = UP-US) = -0.10 to 0.15 V), covering a physiological range of cell membrane potential. [START_REF] Kandel | Principles of Neural Science[END_REF] I-V measurements of half anti-Rhodopsin antibody alone, and hOR1A1 immobilized through half anti-Rhodopsin antibody are shown in Figure 2A. The presence of hOR1A1 caused a pronounced reduction of conductance. We observed a linear I-V response in all cases; thus, a single conductance value (fitted I-V slope) can be assigned to each I-V trace. Figure 2B shows measured distributions of the conductance values obtained for anti-Rhodopsin antibody, and hOR1A1 immobilized through half anti-Rhodopsin antibody with and without the presence of the OR cognate ligand dihydrojasmone at 30 µM concentration. The presence of the ligand caused a significant increase of conductance in hOR1A1.

To exclude contributions from the media in the gap between the probe and the sample, we conducted static break junction experiments. [START_REF] Artés | Transistor-like Behavior of Single Metalloprotein Junctions[END_REF][START_REF] Ruiz | Bioengineering a Single-Protein Junction[END_REF] With the feedback loop transiently disconnected, we recorded current vs. time (I-t) at a constant bias of 0.2 V (US = 0.25 V and UP = 0.45V) and at a current set point of 0.3 nA. Spontaneous contact between the probe and the sample results in jumps (blinks) of the current that last as much as the contact does (Figure S2A). They represent solely the net conductance flowing through the receptor, thus without the contribution of the media. By collecting the blinks with subtracted current baseline and setting them to a common time origin, 2D-blinking maps were built. These maps show the variability in conductance that is introduced by the geometry of the contact (Figure S2B). Taking the most probable conductance peak, we observed that binding to dihydrojasmone caused a conductance increase (Figure S2C) as in I-V non-contact measurements. In general, conductance values obtained in I-t measurements are lower, as expected from the more resistive environment of the protein.

Treatment with increasing concentrations of dihydrojasmone (0-30 µM) lead to a growth of conductance (Figure 2C). Fitting the data using Hill sigmoidal equation with variable slope (i.e., 4 parameters logistical) (Figure S3) yielded an EC50 of 11.2 µM, in agreement with previous results, [START_REF] Belloir | Biophysical and Functional Characterization of the Human Olfactory Receptor OR1A1 Expressed in a Mammalian Inducible Cell Line[END_REF] and thereby validating our EC-STM set-up for studying hOR1A1-dihydrojasmone binding. A Hill coefficient of 9.9 was obtained from the fit indicating a non-linear relation between de EC-STM response and odorant concentration. Supralinear relation (i.e., Hill coefficient > 1) has been previously reported for ORs and attributed to repeated binding. [START_REF] Bhandawat | Elementary Response of Olfactory Receptor Neurons to Odorants[END_REF][START_REF] Sanmartí-Espinal | Quantification of Interacting Cognate Odorants with Olfactory Receptors in Nanovesicles[END_REF] We also conducted current-distance (I-z) measurements on hOR1A1 in the absence/presence of dihydrojasmone at the same concentrations used in I-V experiments (Figure S4). I-z measurements were conducted as previously described, [START_REF] Lagunas | Long Distance Electron Transfer through the Aqueous Solution between Redox Partner Proteins[END_REF][START_REF] Artés | Direct Measurement of Electron Transfer Distance Decay Constants of Single Redox Proteins by Electrochemical Tunneling Spectroscopy[END_REF] at a current set point of 0.4 nA, and at a constant bias of 0.2 V (US = 0.25 V and UP = 0.45V). Up to 100 I-z curves were recorded per sample and distance decay factors (β) were quantified from individual semi-logarithmic I-z plots. We observed that β decreases with the increasing dihydrojasmone concentration, suggesting that dihydrojasmone increases the spatial span of the hOR1A1 currents.

Besides the conductance increase, I-V curves showed a shift of around 20 mV in the open-circuit voltage (VOC) towards lower potentials (absolute values) in the presence of dihydrojasmone (Figure 3A), which is not observed in the absence of the ligand (Figure S5A). The VOC value obtained from the I-V curves for the hOR1A1 was 33.8 ± 0.2 mV, which was similar to that obtained for the half anti-Rhodopsin antibody alone (VOC = 31.5 ± 0.2 mV) and significantly higher (P < 0.0001) than that obtained for hOR1A1 with dihydrojasmone at 30 µM concentration (VOC = 21.8 ± 0.3 mV; Figure 3B). This indicates that ligand binding produces a decrease in the receptor potential. Treatment with increasing concentrations of dihydrojasmone (0-30 µM) caused VOC to decrease exponentially with a 1/e of 12.8 µM (Figure 3C andS5D).

The existence of a Voc different of zero in the experimental results suggests that hOR1A1 acts as an electrical first order system to a voltage ramp. Reported EIS measurements in bulk showed that ligand binding in ORs can be monitored following the variation of the impedance spectra. Nyquist plots were fitted using a modified Randles equivalent circuit, in which the response of the receptor was essentially described by the impedance of the Rp-CPE parallel circuit. [START_REF] Alfinito | Olfactory Receptor-Based Smell Nanobiosensors: An Overview of Theoretical and Experimental Results[END_REF] In agreement with that, we modeled the electrical behavior of hOR1A1 as an RC parallel circuit (Figure 4A), and empiric impedance values were directly determined, resulting in R1 = 583 MΩ, C1 = 0.10 nF (SI). The interaction of OR with dihydrojasmone produces a decrease of around an 11% of the resistive part, and around a 30% of the capacitive part of the electrical equivalent circuit (from 683 MΩ to 517 MΩ, and from 0.10 nF to 0.07 nF, respectively. See SI).

In conclusion, we have measured the nanoscale conductance of hOR1A1 in nearly physiological conditions and we have found that dihydrojasmone binding to the receptor causes an increase in conductance that can be measured from both I-V and I-t recordings, and which is dose dependent. The EC50 of 11.2 µM determined is within the range of EC50 values previously reported for hOR1A1, [START_REF] Belloir | Biophysical and Functional Characterization of the Human Olfactory Receptor OR1A1 Expressed in a Mammalian Inducible Cell Line[END_REF] thus indicating that changes in conductance are correlated with ligand binding. Strikingly, dihydrojasmone binding produces a shift in VOC towards lower potentials. This indicates that ligand binding induces a depolarization of the receptor that causes a decrease in the impedance, in good agreement with the microscopic model proposed by Alfinito et al. for ORs, which is based on the outcome of bulk experiments. [START_REF] Hou | A novel detection strategy for odorant molecules based on controlled bioengineering of rat olfactory receptor I7[END_REF][START_REF] Alfinito | Human Olfactory Receptor 17-40 as an Active part of a Nanobiosensor: a Microscopic Investigation of its Electrical Properties[END_REF] Dihydrojasmone binding causes a conformational change in the receptor that may alter charge distribution facilitating charge transfer, as demonstrated by Iz measurements, and thereby increasing conductance (Figure 4B) and spatially extending (reducing) the β decay rate. The dependence of both conductance and VOC with ligand concentration agrees with the proposed mechanism of activation of GPCRs, for which a complex series of conformationally stable intermediates have been described to evolve towards a more active state with the increasing ligand concentration. [START_REF] Kobilka | Conformational Complexity of G-Protein-coupled Receptors[END_REF] Finally, the simultaneous measurement of RC equivalent by means of the VOC potential allows increasing the electrical sensitivity at single receptor level for biosensing applications. The electric model and fit predicts that faster I-V ramps would provide even larger VOC shifts and thus higher ligand sensitivity both for fundamental and sensing applications.

Overall, EC-STM proved as a reliable technique to study the nanoscopic electrical properties of hOR1A1, which were measured with unprecedented control over the receptor orientation and in nearly physiological conditions. The changes induced by ligand binding in conductance and specially in VOC values, pave the way towards the development of better biohybrid odorant sensors.

ASSOCIATED CONTENT

Supporting Information

The Supporting Information is available free of charge on the ACS Publications website. Materials, sample preparation, details on EC-STM measurements, statistics, atomic force microscopy measurements, conductance fitting, I-z and I-t measurements, and VOC analysis (PDF)

  Lagunas*, Christine Belloir, Loïc Briand, Pau Gorostiza, Josep Samitier* Anna Lagunas -Biomedical Research Networking Center in Bioengineering, Biomaterials, and Nanomedicine (CIBER-BBN), 28029 Madrid, Spain; Institute for Bioengineering of Catalonia (IBEC), Barcelona Institute of Science and Technology, 08028 Barcelona, Spain. Christine Belloir -Centre des Sciences du Goût et de l'Alimentation, AgroSup Dijon, CNRS, INRAE, Université Bourgogne Franche-Comté, Dijon, France. Loïc Briand -Centre des Sciences du Goût et de l'Alimentation, AgroSup Dijon, CNRS, INRAE, Université Bourgogne Franche-Comté, Dijon, France. Pau Gorostiza -Institute for Bioengineering of Catalonia (IBEC), Barcelona Institute of Science and Technology, 08028 Barcelona, Spain; Biomedical Research Networking Center in Bioengineering, Biomaterials, and Nanomedicine (CIBER-BBN), 28029 Madrid, Spain; Catalan Institution for Research and Advanced Studies (ICREA), Barcelona 08011, Spain. Josep Samitier -Institute for Bioengineering of Catalonia (IBEC), Barcelona Institute of Science and Technology, 08028 Barcelona, Spain; Biomedical Research Networking Center in Bioengineering, Biomaterials, and Nanomedicine (CIBER-BBN), 28029 Madrid, Spain; Department of Electronics and Biomedical Engineering, University of Barcelona (UB), 08028 Barcelona, Spain.

Figure 2 .

 2 Figure 2. (A) Ensemble of EC-STM current-bias voltage (I-V) curves obtained for half anti-Rhodopsin (grey) and hOR1A1 (blue), showing the decrease in conductance upon receptor immobilization. (B) Conductance histograms obtained from the linear fitting of individual I-V curves for half anti-Rhodopsin (grey), hOR1A1 before (blue) and after incubation with 30 µM dihydrojasmone (red). (C) Plot of the conductance variation in hOR1A1with the increasing dihydrojasmone concentration. Dashed black lines are an eye-guide. Values are the mean ± SE. n = 150. * P < 0.03, ** P < 0.01, **** P < 0.0001. Experiments were conducted at a constant sample potential (US) of 250 mV, current set point = 0.4 nA, 50 mM sodium phosphate buffer, pH 7.4.

Figure 1 .

 1 Figure 1. (A) Binding activity of purified hOR1A1 using intrinsic tryptophan fluorescence. Dose-response relationship of hOR1A1 fluorescence (λexc = 280 nm, λexc = 340 nm) following dihydrojasmone application. The data were fitted with sigmoid dose-response curves and Kd value obtained (mean ± s.d, N = 3). (B) Schematic representation of the experimental setup. Due to lack of any crystalline structure for hOR1A1, sequence homology with Phyre2 has been used for the receptor representation. Abbreviations: WE, working electrode; RE, reference electrode; CE, counter electrode.

Figure 3 .

 3 Figure 3. (A) Superimposition of I-V curves from hOR1A1 before (blue) and after (red) incubation with 30 µM dihydrojasmone ligand. Dihydrojasmone causes a shift towards a lower absolute value of the open-circuit voltage (VOC), inset. (B) VOC histograms obtained from the linear fitting of individual I-V curves in A. (C) Plot of the VOC variation (absolute values) in hOR1A1with the increasing dihydrojasmone concentration. Dashed black lines are an eye-guide. Values are the mean ± SE. n = 150. * P < 0.03, ** P < 0.01, **** P < 0.0001. Experiments were conducted at a constant sample potential (US) of 250 mV, current set point = 0.4 nA, 50 mM sodium phosphate buffer, pH 7.4.

Figure 4 .

 4 Figure 4. (A) Electrical equivalent model proposed: a combination of RC parallel circuits for the OR (R1C1) and the half antibody (R2C2). The values of R1 and C1 are regulated by the ligand concentration. (B) Dihydrojasmone binding to hOR1A1 may cause charge reorganization facilitating charge transfer and increasing conductance.

ACKNOWLEDGMENT

This work was supported by the Biomedical Research Networking Center (CIBER), Spain. CIBER is an initiative funded by the VI National R&D&i Plan 2008-2011, Iniciativa Ingenio 2010, Consolider Program, CIBER Actions, and the Instituto de Salud Carlos III (RD16/0006/0012; RD16/0011/0022), with the support of the European Regional Development Fund (ERDF). This work was funded by the CERCA Program and by the Commission for Universities and Research of the Department of Innovation, Universities, and Enterprise of the Generalitat de Catalunya (2017-SGR-1079, 2017-SGR-1442 and 2017-SGR-00465). This work was supported by grants from the Conseil Régional Bourgogne Franche-Comté (PARI grant) and the FEDER (European Funding for Regional Economical Development). This research also received funding from the European Union Research and Innovation Programme Horizon 2020 (Human Brain Project SGA3 No. 945539), DEEPER (ICT-36-2020-101016787), Agency for Management of University and Research Grants of the Government of Catalonia (CERCA Programme; 2017-SGR-1442 project; AGAUR, Clúster Emergent del Cervell Humà), and Ministry of Economy and Competitiveness (Grant PID2019-111493RB-I00). The authors also acknowledge Dr. Albert Cortijos Aragonés for help with I-t measurements and data analysis.

AUTHOR INFORMATION

Notes

The authors declare no competing financial interests.


Jean-François Olivieri

Prof Damien J T Laage
email: damien.laage@ens.psl.eu

James T Hynes
email: casey.hynes@ens.psl.eu

Dr J.-F Olivieri

PASTEUR J T Hynes

A Model Electron Transfer Reaction in Confined Aqueous Solution

Keywords: Electron transfer reaction, Confinement, Dielectric constant, Solvent reorganization energy, Aqueous solution

Liquid water confined within nanometer-sized channels exhibits a strongly reduced local dielectric constant perpendicularly to the wall, especially at the interface, and this has been suggested to induce faster electron transfer kinetics at the interface than in the bulk. We study a model electron transfer reaction in aqueous solution confined between graphene sheets with classical molecular dynamics. We show that the solvent reorganization energy is reduced at the interface compared to the bulk, which explains the larger rate constant. However, this facilitated solvent reorganization is due to the partial desolvation by the graphene sheet of the ions involved in the electron transfer and not to a local dielectric constant reduction effect.

Introduction

The impact of confinement on structural and dynamical properties of liquids is of interest in a broad variety of contexts, including, e.g. transport within nanofluidic devices to catalysis of chemical reactions in nanoreactors [1][2][3][4][5] . In one of these, recent experimental [6,7] and simulation [8][9][10][11][12] studies have suggested that the dielectric constant of water is dramatically reduced in nanometer-size channels along the direction perpendicular to the channel wall, and exhibits very strong anisotropy with respect to the confining media. Such effects could have considerable importance for charge transfer reactions: such reduced permittivity could facilitate transfer by reducing solvent reorganization penalties which slow such reactions. Indeed, it has been suggested that very rapid electron transfer (ET) reaction rates at interfaces determined in recent experimental measurements could be due to the reduced interfacial dielectric constant. [13] In the present work, we explore such a possibility for key aspects of a model ET reaction 27 in water bounded by (here uncharged) graphene sheets, whose separation is fixed, a system of considerable interest in a wide variety of arenas. [14][15][16][17] Our study employs both molecular dynamics simulations and the theoretical results of previous efforts [12,18,19] on water confined between such sheets, in particular a recent investigation [12] of water's dielectric constant behavior in such an environment. Many prior simulation studies (see e.g. [20][21][22][23][24][25][26][27][28][29][30][31][32] ) have focused on ET between a redox species in solution and a metallic electrode. These have identified that for these systems, the major effect on ET reaction is due to the interface's metallic nature. The image charge effect due to electronic polarization in the metal reduces the solvent reorganization energy λ [28,32] , accelerating the ET reaction, an effect which can extend several nm from the interface [28] . In addition, λ may be influenced by modified interfacial solvation shell structure. [21,26,30,33] Here we are instead concerned with an aqueous solution model ET reaction confined between graphene surfaces, which do not involve such significant image charge effects; this allows a focus of attention on interfacial region (and beyond) solvation effects, including those related to water's dielectric constant behavior. In this first effort, we focus on an ET reaction not involving charge transfer to or from the graphene, but rather examine an ET for a solute within the water solvent. Further, we do not address all aspects of the ET reaction rate, concentrating instead on the reorganization energy λ which has an important impact on that rate [34][35][START_REF] Newton | Electron transfer in chemistry[END_REF][START_REF] Savéant | Elements of molecular and biomolecular electrochemistry[END_REF] . The model employed [START_REF] Zichi | [END_REF] is quite simple and, as discussed within, has a number of simplifying features which allow this focus on solvation effects unencumbered by other important, but also obscuring, features influencing the reaction. In this fashion we are able to elucidate the role of several featuresincluding the water dielectric constant -governing the behavior of λ both close to and away from the graphene surfaces. In particular, we show that the reduced solvent reorganization energy for ET at the graphene interface is not due to the reduced dielectric constant component perpendicular to the wall. Instead, we show that the small reorganization energy results from the partial desolvation induced by the wall for solutes lying in the first layer. The outline of the remainder of this paper is the following. We first describe the model electron transfer reaction and introduce the key activation free energy and solvent reorganization energy quantities, followed by a presentation of our simulation methodology. We then analyze the free energy curves and reorganization energies obtained from our simulations in different geometries for the solute within the confined water slab, with a special focus on the role of the effective dielectric constant. Finally, we offer some concluding remarks.

Model electron transfer reaction

The model ET reaction system in graphene plate-confined water selected for this molecular dynamics study is closely related to that introduced by Zichi et al [START_REF] Zichi | [END_REF] . It involves two atomic charge centers A and B separated by a variable distance RAB, the reaction consisting of the single electron transfer

A -1/2 B +1/2 ® A +1/2 B -1/2
(1) This ET has the advantage that it is thermodynamically symmetric, i.e., the reaction free energy ΔGrxn = 0, in the bulk. This both avoids any influence of the reaction thermodynamics on the ET barrier in the bulk and, most especially, allows a clear exposure of any reaction asymmetry which may be introduced in the interfacial region by the presence of the graphene surfaces. In contrast to Ref. [START_REF] Zichi | [END_REF] , we do not address in our study the electronic coupling Vel between the diabatic reactant and product valence bond states, which would (strongly) vary with the A-B separation RAB and importantly impact the ET reaction rate: lowering the barrier in the adiabatic regime (smaller RAB, larger Vel) and providing a rate prefactor in the nonadiabatic tunneling regime (larger RAB, smaller Vel). Omitting these impacts allows an exclusive focus on the solvent reorganization energy λ, which importantly contributes to the rate whatever the value of Vel (cf eq. 2). On the other hand, in order to clarify the significance of key quantities that we calculate, it is useful to quote the well-known Marcus equation for a nonadiabatic, small Vel ET reaction barrier ΔG ‡ involving the already mentioned reorganization energy and the reaction free energy [34][35][START_REF] Newton | Electron transfer in chemistry[END_REF][START_REF] Savéant | Elements of molecular and biomolecular electrochemistry[END_REF] ΔG ‡ = (λ + ΔGrxn) 2 / 4 λ (2) which shows that the reaction barrier is λ/4 for the thermodynamically symmetric reaction. For a symmetric adiabatic reaction, this barrier would be reduced by Vel. [35,[START_REF] Newton | Electron transfer in chemistry[END_REF]39] Equation . 2 relies on a certain harmonic approximation discussed in a subsequent section, where relevant free energy curves are found via simulation rather than harmonic approximation. Also useful for perspective in our efforts is (in a form appropriate for our model) the Marcus model dielectric continuum approximation for λ [34][35][START_REF] Newton | Electron transfer in chemistry[END_REF] λ

(R) = Δq 2 (1-1/ε) (1/2rA + 1/2rB -1/RAB) (3)
where Δq = 1 e is the transferred charge, ε is the solvent dielectric constant (with electronic polarization neglected), rA,B are the radii of the A and B atoms, and RAB is the A-B separation. Equation 3 reflects the generally expected trends in a bulk solution symmetric reaction: λ -and thus the ET barrier -decreases as ε decreases towards 1, and increases to an asymptotic value as the A-B separation RAB increases. (A feature that will be important for our analysis is that this very large RAB asymptotic limit shares the same dependence on ε and on rA,B as the familiar Born model [40] solvation free energies ΔGsolv Born for A and B: e.g., eq. 3 reduces to the sum of these free energies, scaled by (Δq/qA,B) 2 ; in addition, eq. 3 assumes an isotropic environment while this is not the case at the interface, and we will return to this point in the following) The eq. 1 model's van der Waals radii are identical for A and B, but the more relevant radii for eq. 3 are the cavity radii, which can be approximated by the location of the first peak in the solutesolvent radial distribution function. This peak occurs at a shorter distance for the anion than for the cation (see Fig. 1f), a wellknown asymmetry for water solvent [41][42][43] . Thus, in the bulk the ET reaction is symmetric, although the solvation free energies for the reaction complex's two charged moieties differ. In the simulations, several different geometries will be considered (Scheme 1): i) A is within the 1 st layer next to graphene and B is located such that the AB vector is normal (perpendicular) to the wall; ii) Both A and B are located within the 1 st layer next to graphene, such that the AB vector is parallel to the wall; for both geometries, the A-B separation RAB is varied, with special attention paid to the reference situations (iii) and (iv), in which the perpendicular and parallel geometries are also examined at the midpoint of the water slab. Finally, the A-B system is separately considered in the bulk, for reference.

Methodology

Classical molecular dynamics simulations are performed with the model ET reactive solute (described above) immersed in SPC/E water within a 6 nm-thick periodic slab between graphene sheets, as well as in a bulk SPC/E water reference system. The two 240atom rigid graphene sheets are constructed from an ideal hexagonal lattice with a 1.42 Å distance between neighboring carbon atoms. There is zero charge placed on the graphene carbon atoms (see Refs. [18,19] for the distinct situation of graphene electrodes held at constant potential and the impact on interfacial water layer structure and dynamics). The simulation box dimensions are 24.58x25.56x62 Å 3 and the number of water molecules ( 1150) is determined such that the water has the correct bulk density in the middle of the slab. The Lennard-Jones parameters for sites A and B are σA,B = 2.5 Å and εA,B = 100 kBT (additional results for the larger solute parameter σA,B = 3.5 Å are reported in the SI). Graphene carbon Lennard-Jones parameters are taken from ref. [44] where they were optimized to reproduce DFT-based molecular dynamics simulation results (see SI for details). Water molecules are held rigid with the SHAKE algorithm. The Yeh-Berkowitz slab correction [45] to Ewald sums for 2D systems is employed. Scheme 1. Different geometries for the A-B pair within the water slab between the graphene plates.

For each system, five independent initial configurations are prepared. Molecular dynamics simulations are run with LAMMPS [46] . After 1 ns equilibration at 300 K, 1 ns production runs are propagated with a 2 fs timestep and the velocity Verlet integrator in the canonical ensemble, using the velocity rescale thermostat [47] with a 2 ps relaxation time. Configurations are output with a 50 fs interval. Solute-solvent interaction energies are calculated with both reactant and product charge distributions at every output step. This procedure is repeated for each geometry shown in Scheme 1 and for a series of A-B separations by freezing the positions of A and B. The bulk reference system trajectories were produced following the same procedure, with 40 independent starting configurations, using a 30x30x30 Å 3 box containing 898 water molecules and a pair of A and B solutes whose positions are not constrained. The free energy profiles of the reactant (R) and product (P) diabatic states are calculated from a series of independent

Perpendicular Parallel Interface Middle i) ii) iii) iv) R AB B A R AB A B R AB A B R AB B A
trajectories propagated with the mapping potential Vη = (1-η) VR + η VP, with η ranging from 0 (pure state R) to 1 (pure state P). The mapping potential is used in this work to generate solvent configurations equilibrated to different electronic structures of the AB pair. The free energy of each diabatic state is reconstructed using the weighted histogram analysis method (wham) [48] . While wham is usually employed with harmonic biasing potentials, it can be easily generalized to the present mapping potential simulations by defining the biasing potential as wη R (ΔE) = Vη(ΔE) -VR(ΔE) = -η ΔE, and wη P (ΔE) = (1-η) ΔE. (Here ΔE is a solvent coordinate defined at the beginning of the following section.) The diabatic free energy profiles are determined using 11 independent 1 ns-long windows regularly spaced along η. The solvation free energy profile of solutes A and B at varied separations from the wall is calculated using a series of independent trajectories where a single solute is present in the box and its distance z to the wall is harmonically constrained by k (z-z0) 2 with k = 50 kBT/Å 2 and z0 ranges from 1.5 Å to 30 Å in 0.125 Å increments. Five independent 200 ps trajectories are used for each z0 distance and the free energy profile is reconstructed with the wham [48] algorithm.

ET Free Energy Curves

Reaction Coordinate

We need to specify a reaction coordinate in order to generate diabatic free energy curves for the ET reaction to determine the reorganization energy λ as well as the reaction free energy ΔGrxn. Here we select a reaction coordinate ΔE is designed to allow the sampling of different, usually nonequilibrium, arrangements of the reaction solute's environment. ΔE is defined as a certain difference in the potential energy of interaction of the different charge states of the reaction solute complex A-B with the water molecules and the graphene plates (eq. 4). We now anticipate that in our model -as to be explained presently -the potential energies in eq. 4 are simply the Coulomb potential interaction between the charges in the reaction solute sites and those of the water solvent molecules. Then ΔE(x) = VR(x) -VP(x) (4) is the difference of the Coulomb interaction energies of the solvent water molecules with the different varied charge distributions of the reactant R complex VR(x) and of the product P complex VP(x), where x designates the solute and solvent positions. The detailed description of the calculation of the free energy curves in this coordinate is given in the Methodology section. We now return to the postponed issue of eq. 4's important simplification that only Coulomb reaction complex-water molecule interactions are involved. This results from our reaction model and our restriction to uncharged graphene plates: the non-Coulombic interactions of the model's A and B units with the solvent waters and the graphene plates remain the same independently of the charges on them, and thus completely vanish in the difference eq. 4.

Free energy curves

Figure 1 shows the diabatic free energy curves calculated from our simulations in a series of geometries for the AB pair (see Scheme 1) and for RAB distances of 5 Å (which approximately corresponds to a solvent separated ion pair arrangement, see Fig. 1f) and 10 Å (beyond which the electronic coupling exponential decay with the separation will make the ET reaction very slow). In the following, we analyze these free energy curves and their harmonic approximations, and successively discuss the reaction free energy, reaction free energy barrier and solvent reorganization energy. ). Solid lines show the profiles obtained from a series of simulations biasing the solvent coordinate between its two equilibrium values for RAB=5 Å (blue) and RAB=10 Å (red). Dashes are the parabolic profiles obtained within the harmonic approximation from the equilibrium solvent coordinates (see e.g. [35,49] ). f) Radial distribution function between anionic solute A and water oxygen atoms (blue), cationic solute B and water oxygen atoms (red), and between A and B (green) in a bulk aqueous solution.

Equilibrium Reaction Free Energy ΔGrxn

A first aspect of the preceding free energy curves to discuss is the equilibrium reaction free energy ΔGrxn. In the bulk, in the middle of the slab and in the parallel arrangement at the interface, ΔGrxn=0 by symmetry. This contrasts with the situation at the interface with the perpendicular geometry (recall that A remains in the interfacial region while B moves away). The free energy profiles in Fig. 1a show that the ET reaction is slightly endergonic: ΔGrxn = 0.55 kcal/mol when RAB=5 Å, and ΔGrxn = 0.45 kcal/mol when RAB=10 Å (note that the harmonic approximation incorrectly predicts a strongly exergonic ET reaction for both RAB distances; we return to that overestimate in our discussion of Fig. 4).

The key to understanding this free energy cost is the feature that the eq. 1 ET removes a negative charge from the surfaceadjacent A, replacing it with a positive charge. This strongly suggests that ΔGrxn > 0 results from the greater stability of the anion in the interfacial region compared to that of the cation. This is confirmed in Fig. 2 which displays the solvation free energy profiles for the separate anion and cation as a function of distance z from the graphene plate. By the same token, if the orientation of the reactant pair were reversed, with B being the surface-adjacent member of the complex, the interfacial ET reaction becomes exergonic ΔGrxn < 0. The propensity of ions to lie in the bulk or at a hydrophobic aqueous interface (including, e.g., the air/water interface) has been shown [50] to be determined by the strength of the ion-water interaction (and by ion polarizability, not included in our simulations). The calculated solvent coordination numbers around the anion and cation solutes for different geometries and RAB separations (Fig. 3) provides a molecular understanding of the reason why the anion is more stable than the cation in the first layer next to the apolar graphene wall. The anionic solute strongly interacts with the water hydrogen atoms, which do not possess any Lennard-Jones sphere and can approach close to the anion. This explains the smaller cavity radius for the anion than for the cation (Fig 1 .f) and the smaller coordination number for the anion than for the cation (Fig. 3). When the solutes are moved from the middle of the slab to the interface, Fig. 3 shows that the decrease in coordination number is much more pronounced for the larger cation than for the smaller anion. The graphene wall-induced partial desolvation thus leads to a larger free energy cost for the cation than for the anion. The difference in solvation free energies for the cation and anion at the graphene interface (≃0.6 kcal/mol) is in quantitative agreement with the ΔGrxn reaction free energy asymmetry.

ET Solvent Reorganization Energy and Reaction Barrier

In order to determine the impact of the interface and of the AB pair orientation on the ET rate constant, we now turn to the solvent reorganization energy λ and the closely related reaction barrier. We recall from the Introduction and section "Model electron transfer reaction" that λ is key for the ET reaction rate via its barrier, which increases with increasing λ. We first determine λ in each diabatic state from our simulated free energy profiles shown in Fig. 1 as λi = Gi(ΔEj eq ) -Gi(ΔEi eq ) ( 5) where i,j=R,P, and the ET free energy barrier as ΔG ‡ = GR(0) -GR(ΔER eq ) ( 6) We then compare these values with the harmonic approximation (sometimes referred to as "linear response") [49,51] λharm = (ΔEP eq -ΔER eq ) / 2 (7) and eq 2 for the free energy barrier in the harmonic approximation. The results in Fig 4 show that λ and ΔG ‡ are very sensitive to the local environment (interface, middle of the slab, bulk) and to the AB pair orientation: λ and ΔG ‡ are smaller at the interface than in the middle of the slab, and this reduction is even more pronounced for the parallel geometry than for the perpendicular one. The fairly large solvent reorganization energies in Fig. 4a arise from the model system's small solute radii and the resulting strong solute-solvent interaction; while these values are not meant to be realistic, this helps the simple harmonic approximation, which will facilitate our analysis (as shown in the SI, our conclusions are unchanged when a solute pair with weaker solute-solvent interaction is considered). In addition, the large reorganization energies λ imply that when the A-B complex is within the interfacial region, the ΔGrxn departure from a symmetric reaction (comparable to the thermal energy kBT and larger than this, though overestimated, in the harmonic approximation, see Fig. 1) has a small impact on the simulated reaction barrier height (Fig. 4b). A related consequence is that the Marcus inverse regime is not observed for this system. Figure 4a shows that the harmonic approximation provides a good description of the changes in λ between the interfacial, middle and bulk environments, and between the parallel and perpendicular orientations. One exception where the harmonic approximation deviates more markedly from the simulation results is the perpendicular geometry at the interface, due to the substantial difference between the λR and λP reorganization energies around the R and P states. As shown by the coordination numbers in Fig. 3, this is caused by the contrast between the solvation structures in the first interfacial layer and in the further layers. The harmonic approximation also provides a good description of the reaction free energy barrier ΔG ‡ (Fig. 4b) and its changes with the different environments and geometries. Its main limitation is its systematic but limited underestimation of the reaction's free energy, due to its underestimation of the solvent reorganization energy in the reactant state (Fig. 4a).

For simplicity, we adopt the harmonic approximation in the following, and we take advantage of its reduced computational cost compared to the full free energy profile calculation to examine how λ changes over a broader range of locations for the reaction A-B complex. Figure 5 displays the harmonic approximation results eq 7 versus the separation RAB of the A-B pair for the pairs' parallel and perpendicular orientation with respect to the graphene surface, successively for the complex with its A member always fixed in the first layer at the graphene interface, then for A in the second layer, for the complex centered in the middle of the two graphene plates, and finally for the complex in the bulk. In all cases, λ increases with increasing RAB. At smaller RAB, the hydration shells of the anionic A and cationic B complex members interfere with each other to a degree. As RAB increases, λ will increase, asymptotically approaching the higher cost of the reorganization of the separate, full hydration shells of the anionic and cationic components of the complex [52] . On the way to this limit, in the middle of the 6 nm-thick slab, the orientation of the reaction complex's orientation matters less for λ than at the interface; a residual anisotropy is nonetheless present, consistent with the long-range (approximately 3 nm) decay length of the interface's effect on the dielectric constant along the direction perpendicular to the wall that was determined in our recent study [12] . On the other hand, in the interfacial geometries where the reorganization cost is lower, the parallel geometry λ is noticeably lower than that for the perpendicular case. Finally, in the perpendicular geometry, moving A from the first layer to the second layer next to the graphene wall leads to λ values which are very close to those in the middle of the slab. We will return to these features presently.

Dielectric Constant and Reorganization Energy

We now determine the molecular origin of the λ reduction at the interface, and examine the suggestion [13] that it is due to the reduced local dielectric constant perpendicularly to the wall. Simulations have shown that the local dielectric constant at an interface is highly anisotropic: it is reduced relative to the bulk along the direction perpendicular to the wall, and typically enhanced along the parallel direction. This anisotropy is included in our simulations and is for example visible in the differences between the free energy profiles obtained with parallel and perpendicular orientations in Fig. 1. To analyze the changes in the solvent reorganization energy at the interface, we employ eq. 3.

The latter was derived with a number of approximations, including a continuum solvent description and an isotropic environment; regarding the latter, its tensorial extension could be envisaged.

Our focus is however on key aspects of the ET reaction rate. We therefore use eq. 3 and probe the system's anisotropy by examining λ along the directions parallel and perpendicular to the wall and determining the impact of the dielectric constant anisotropy on the solvent reorganization energy. To this end, we consider the different solute pair orientations described in Scheme 1, and use the simple isotropic continuum model eq. 3 to probe the effective dielectric constant for each orientation. We note that, remarkably, in all cases in Fig. 5, the A-B separation dependence of λ found in the molecular simulations is very close to the 1/RAB dependence which follows from the dielectric continuum model eq. 3 (see also Fig. 6 below). In eq 3, the slope of λ with 1/RAB only depends on 1/ε (and on Δq which is fixed in our model eq. 1). On the other hand, the intercept of λ along 1/RAB depends both on 1/ε and on the cavity radii rA;B whose values cannot be unambiguously determined from the molecular simulations. with the perpendicular (cyan) and parallel (orange) geometries (see Scheme 1). b) Reaction activation free energy from the free energy profiles as a function of the harmonic activation free energy (eq 2), with the same symbols as in panel a.

Figure 5. Reorganization energy (eq 7) as a function of the RAB distance plotted with a reciprocal scale for the horizontal axis, respectively in the perpendicular (black) and parallel (red) geometries at the interface, in the perpendicular (green) and parallel (blue) geometries in the middle of the slab, in the perpendicular geometry with A in the second layer next to the graphene sheet (pink) and in the bulk (orange) (see Scheme 1) (the pink and green symbols strongly overlap). The chosen 3-10 Å range for RAB extends from typical A-B distances for a contact ion pair (Fig. 1f) to long-range ET (limited in our simulations by the periodic simulation box size)

Given this slope dependence, one could try to determine an effective ε from eq. 3 and the slope of λ versus 1/RAB. Unfortunately, due to the inverse proportionality in the dielectric constant, we find that the resulting uncertainty in ε is too large to obtain meaningful ε values. We therefore adopt a different strategy. To analyze the dielectric constant aspects of the reorganization energy, we focus on the interfacial case, with the reaction complex member A always in the interfacial layer, with member B at various separations RAB from it and with the AB direction either perpendicular or parallel to the graphene sheet. In particular, we compare the simulation results for λ with the dielectric continuum formula eq. 3, with three different choices for the description of the dielectric constant ε: i) the bulk SPC/E dielectric constant ε = 71 [53] , ii) the local interfacial dielectric constant ε ≃ 2 that was suggested [7,10] from a 2-state capacitor series interpretation of experimental and simulation results, assuming an extremely restricted first layer and a bulklike behavior in further layers [START_REF]While ε =2 is the continuum value for a solely electronically polarizable liquid, our molecular dynamics water solvent is not[END_REF] , and iii) the local dielectric constant profile that we recently determined [12] from a Kirkwood factor analysis of molecular dynamics simulations. This last picture is sufficiently novel to require a brief description here. The dielectric constant along the direction perpendicular to the wall is reduced at the graphene interface due to the anisotropic dipolar field and the exclusion by the wall of water molecules that would bring a positive dipole-dipole correlation in the bulk. The dipoledipole interaction's long-range character induces a slow decay of the dielectric constant reduction on a scale of several nanometers from the wall. Due to this effect's strong anisotropy, the dielectric constant along the direction parallel to the wall is much less affected by the dipolar field anisotropy and excluded-volume combination; simulations showed that it is actually locally enhanced compared to the bulk [9,10,12,[START_REF] Motevaselian | [END_REF] (ε = 175 in the first layer), due to the orientation of the first layer molecules which tend to lie with their dipole tangent to the wall (and only the first layer molecules are affected).

Figure 6a shows, for the perpendicular AB orientation, the 1/RAB dependence of λ both from our simulations and from the dielectric continuum expression eq. 3 with the above-mentioned three different dielectric constant values (the intercept value is numerically adjusted via a linear regression). The ε = 2 case leads to a serious underestimate of the slope and is not compatible with the simulation results. The bulk dielectric constant provides a better description, but the best agreement is reached with the local dielectric constant profile determined from our recent simulations [12] , and whose character was sketched above [START_REF] However | we note that for considerably narrower slabs[END_REF] .

In the opposite, parallel arrangement, Fig. 6b shows again that the dielectric continuum expression with ε = 2 cannot describe the 1/RAB dependence of λ. Both ε = 71 and ε = 175 (cf SI of ref [12] ) provide satisfactory descriptions of the simulation results and the inverse dependence of the slope with the dielectric constant does not allow a more precise determination of its value (However, while the increased ε = 175 value could be consistent with our calculated 1/RAB dependence of λ, it cannot explain the reduction in λ at the interface, since eq. 3 shows that it would instead lead to an increase in λ).

The results in Fig 6 therefore show that the interfacial λ is not consistent with a dramatically reduced dielectric constant at the interface, but that our recent analysis of the local dielectric constant provides an excellent description of the effective dielectric constant entering in the continuum model expression eq. 3.

Since we have concluded that the reduction in λ at the interface is not due to an extremely low dielectric constant, we return to Fig. 5 to determine its molecular origin. Figure 5 shows that while the curves of λ versus 1/RAB exhibit fairly similar slopes across the different environments and orientations, they are vertically shifted, i.e., they have different intercepts (Figures S2-3 show that a solute pair with a weaker solvent interaction yields very similar results). This indicates that the λ reduction on passing from the bulk to the interface is mostly due to a change in the solvation of the AB species. The decrease in coordination numbers shown in Fig. 3 for A and B upon movement from the middle of the slab to the interface reflects a partial desolvation due to the wall, which reduces the energy cost to change that solvation, i.e., the observed solvent reorganization energy decrease at the interface. This conclusion is further supported by the enhanced decrease in λ for the parallel geometry, when both A and B are partly desolvated in the first layer. This stresses the importance of the interfacial solvation shell structure that has been recognized in other contexts [21,26,30,33] . As reviewed e.g. in ref. [START_REF] Benjamin | [END_REF] , this highlights an important limitation of dielectric continuum models which do not correctly describe these hydration structure changes at interfaces. Finally, the considerations of the present work should be relevant for a number of interfacial elementary chemical processes, including not only ET reactions but also proton transfer reactions [58] and solvation dynamics [59] for which reduced solvent polarity at an interface has been shown to have an important impact. Reorganization energy (kcal/mol) Figure 6. a) Reorganization energy (eq 6) in the perpendicular geometry at the interface as a function of the RAB distance plotted with a reciprocal scale for the horizontal axis from our simulations (black symbols) and fits using the continuum equation 3 with ε=71 [53] (green), ε=2 (blue) and the dielectric constant profile ε (z) determined in ref. [12] (red); the ε=2 curve is extended with dashes when B moves beyond the suggested [7] restricted interfacial layer thickness. b) Idem in the parallel geometry, where the red curve is the enhanced ε=175 parallel interfacial dielectric constant determined in ref. [12] , and the green curve ε=71 overlaps with the red curve.

Concluding Remarks

Here we have examined, by simulations and extensive analysis, a model electron transfer reaction within liquid water confined between two graphene sheets. The chosen model system's simplicity allows our focus on solvation and solvent dielectric constant effects for the reaction without the influence of other obscuring features. Our molecular dynamics simulations are analyzed with the Marcus model dielectric continuum approximation for λ and the system's anisotropy is probed by considering parallel and perpendicular orientations of the reactive pair relative to the wall. Our principal conclusion is that a reduction of the reorganization energy and the reaction's activation free energy barrier in the interfacial region is not due to a muchreduced dielectric constant perpendicularly to the graphene sheet, but rather to a wall-induced effect of desolvation of the reaction complex. Future efforts will address improved realistic models of the reaction, assessment of dynamical effects on the rate [START_REF] Zichi | [END_REF] , and most especially, corresponding considerations of electrochemical reactions involving electron transfer to and from an electrode.

Entry for the Table of Contents

We have studied, via simulation and analysis, the reorganization energy and free energy reaction barrier for a model electron transfer reaction occurring at various locations within liquid water confined between two graphene sheets. We find that the acceleration of the reaction in the interfacial region by reduction of both these quantities is due not to significant local dielectric constant reduction, but to desolvation of the reaction complex.

Institute and/or researcher Twitter usernames: @Chimie_ENS @INC_CNRS e -e -

Figure 1 .

 1 Figure 1. Diabatic free energy curves for the AB pair at the interface in the perpendicular (a) and parallel geometries (b), in the middle of the slab in the perpendicular (c) and parallel geometries (d) and in the bulk (e) (see Scheme 1). Solid lines show the profiles obtained from a series of simulations biasing the solvent coordinate between its two equilibrium values for RAB=5 Å (blue) and

Figure 2 .

 2 Figure 2. Solvation free energy profiles for the separate anion (violet) and cation (green) as a function of distance z from the graphene plate.

Figure 3 .

 3 Figure 3. Coordination numbers, defined as the number of water oxygen atoms within the first shell around solute A-B's site A in its anionic reactant (a) and cationic product (b) state as a function of the RAB distance, respectively in the perpendicular (blue) and parallel (red) geometries (see Scheme 1) and where site A lies within the graphene interfacial layer (open symbols) and in the middle of the slab (full symbols).

Figure 4 .

 4 Figure 4. a) Reorganization energy from the free energy profiles (eq 5) as a function of the harmonic reorganization energy (eq 7), in the reactant (full symbols) and in the product (open symbols) diabatic states, at RAB distances of 5 Å (circles) and 10 Å (diamonds), in the bulk (green), at the interface in the perpendicular (blue) and parallel (red) geometries, and in the middle of the slab

Supplementary Information

A Model Electron Transfer Reaction in Confined Aqueous Solution

Jean-François Olivieri, [a] Damien Laage* [a] and James T. Hynes* [a,b] [a] PASTEUR, Department of Chemistry, École Normale Supérieure, PSL University, Sorbonne Université, CNRS 75005 Paris, France E-mail: damien.laage@ens.psl.eu [b] Department of Chemistry, University of Colorado Boulder, CO 80309, USA E-mail: casey.hynes@ens.psl.eu

Simulation Methodology: Combination Rules

For the carbon-oxygen and carbon-solute Lennard-Jones (LJ) interaction potentials, we followed the approach of ref. [1] where carbon-oxygen LJ parameters were optimized to reproduce DFT-based molecular dynamics simulation results. A modified LJ potential was introduced where the attractive term was scaled

and Lorentz-Berthelot mixing rules were used to determine the cross parameters for the interaction between atoms i and j, . !" = ". !! + . "" $/2 and 2 !" = 3 2 !! 2 "" . The optimized values for the graphene carbon atoms were determined in ref. [1] to be εCC=0.3598 kJ/mol and σCC=3.4 Å, and the same tuning parameter λ=0.7 was used for the carbon-oxygen and carbonsolute LJ interactions. Solute-oxygen and solute-water LJ interactions employed the Lorentz-Berthelot mixing rules with no scaling.

Density Profiles

Figure S1 shows the water oxygen and hydrogen number density profiles across the slab, together with a typical configuration of the A-B pair in configuration i (see Scheme 1) with A in the first layer and the AB direction perpendicular to the interface. The relative sizes of the A and B solutes with respect to the hydration layers' thickness are rendered by the colored disks, whose diameters are the Lennard-Jones σA and σB.

Solute Size Effect

We repeated the solvent reorganization energy calculations for larger A and B solutes (σA=σB=3.5 Å), which lead to weaker solute-solvent interaction energies and smaller solvent reorganization energies. Figures S2 andS3 show that the results for these weaker solute-solvent interaction energies exhibit the same behavior for the solvent reorganization energy dependence on the RAB distance as reported in the main text in Figures 5 and6.

Figure S2.

Analog of Figure 5, here for larger A and B solutes (σA=σB=3.5 Å). Reorganization energy (eq 7) as a function of the RAB distance plotted with a reciprocal scale for the horizontal axis, respectively in the perpendicular (black) and parallel (red) geometries at the interface, in the perpendicular (green) and parallel (blue) geometries in the middle of the slab.

Figure S3

. Analog of Figure 6a, here for larger A and B solutes (σA=σB=3.5Å). Reorganization energy (eq 6) in the perpendicular geometry at the interface as a function of the RAB distance plotted with a reciprocal scale for the horizontal axis from our simulations (black symbols) and fits using the continuum equation 3 with ε=71 [2] (green), ε=2 (blue) and the dielectric constant profile ε (z) determined in ref. [3] (red). A fit with the enhanced ε=175 parallel interfacial dielectric constant determined in ref. [3] is indistinguishable from the ε=71 fit (green).

[1] M. K.