Vivien Walter

Céline Ruscher

Olivier Benzerara

Carlos M Marques

Fabrice Thalmann

A machine learning assessment of the two states model for lipid bilayer phase transitions

We have adapted a set of classification algorithms, also known as Machine Learning, to the identification of fluid and gel domains close to the main transition of dipalmitoyl-phosphatidylcholine (DPPC) bilayers. Using atomistic molecular dynamics conformations in the low and high temperature phases as learning sets, the algorithm was trained to categorize individual lipid configurations as fluid or gel, in relation with the usual two-states phenomenological description of the lipid melting transition. We demonstrate that our machine can learn and sort lipids according to their most likely state without prior assumption regarding the nature of the order parameter of the transition. Results from our machine learning approach provides strong support in favor of a two-states model approach of membrane fluidity.

Phospholipid molecules play a major structural role in biological membranes [1,22] where a deep understanding of the physical properties can only be acquired from a detailed knowledge of the lipid assemblies. Thanks to their amphiphilic nature and geometrical characteristics, most phospholipid molecules spontaneously self-assemble in water as bilayers, to eventually form closed vesicles. Lipid vesicles can easily be made, controlled and studied, and have now become standard tools in membrane biophysics studies, referred as model lipid bilayer systems. Early studies on pure phospholipid bilayers indicated that lipids were subject to thermodynamic transitions [3][4][5]22], with in particular a sharp transition associated to a significant change in enthalpy called main, or melting transition. This transition separates a low temperature well-packed assembly from a high temperature expanded, disordered lipid tail organization. This transition is considered to be of first order, although phase coexistence is difficult to establish experimentally [6]. The low temperature phase is commonly referred as the gel state, while the high temperature state is the fluid state. Lipid mixtures also display melting transitions spreading along a finite temperature range, usually accompanied by gel-fluid domain coexistence. It is usually assumed that most biological lipid membranes are found in a fluid state, and many scenarii aiming at explaining the lat- * fabrice.thalmann@ics-cnrs.unistra.fr eral lipid and protein segregation observed in biological membranes involve ordering of the lipid tails.

The consensual description of the single component lipid melting transition assumes that dominant molecular conformations evolve from all-trans extended, well oriented hydrocarbon chain conformations in the low temperature phase, to disordered chains melted by rotation isomerism, as proposed in the earliest theoretical proposals [7][8][9][10]. Lipids in the fluid phase have more configuration entropy and more enthalpy than those in the gel phase, due to lower density and cohesive energy, and higher chain torsion energy. Balance between entropy and enthalpy holds at the coexistence temperature T m . Despite some asymmetry between the two phases, a large number of experimental facts related to melting transition of pure and mixed lipid compositions have been successfully interpreted by means of a phenomenological two-states model, originally proposed by Doniach [4,[11][12][13][14][15][16][17][18]. This model can be expressed as an Ising model, each lipid taking binary discrete values (say s = 0 for gel and s = 1 for fluid ) with neighboring lipids being coupled [4]. In the framework of the two-states model, the Ising variables stand for a coarse-grained description of the lipid tail conformations, assuming that lipids can be classified into two classes, according to their molecular conformations. Within this description, an effective temperature dependent "magnetic field" h(T ) biases the odds in favor of one or the other state, while cooperativity results from nearest neighbor state coupling.

Usually, the determination of the lipid bilayer state relies on a structural scalar order parameter, such as the membrane thickness, given for instance by the head to tail lipid extension or the tail molecular order parameter.

We address in this work the validity of the two-states description for the transition of pure DPPC (1,2-dipalmitoyl-sn-glycero-3phosphatidylcholine) bilayers, using atomistic molecular dynamics (MD) simulations and supervised Machine Learning (ML) classification algorithms. The advantage of using ML classification is to bypass entirely the need of defining a scalar order parameter prior to sorting the lipids into their respective states. Conversely, the ML approach may serve as testing the relevance of a scalar order parameter a posteriori. Machine Learning has already been applied successfully to a number of situations in statistical thermodynamics and phase transitions. For instance, Cubuk et al. used support vector machines to localize plastic flow regions in amorphous structures [19], Carrasquilla and Melko revealed the strong aptitude of neural networks models to recognize various spin ordering regimes in condensed matter systems [20], Le and Tran succeeded in predicting the polymorphism of complex lipid mixtures given a set of structural, chemical and composition parameters, by means of an artificial neural network approach [21]. We show in the present work that Machine Learning classification tools are also extremely powerful for investigating thermal lipid phase transitions. DPPC molecules are among the best known phospholipids [5]. They display experimentally a melting transition at 314 K, which is well reproduced by the CHARMM-36/TIP3P force-field for atomistic simula-tions of lipid bilayers in aqueous solutions [3,22]. As in experiments, the high temperature state is a L α fluid disordered bilayer, and the low temperature state a P β rippled phase (Figure 1), i.e a spatially modulated gel phase with peristaltic variations in the bilayer thickness and normal direction [24,25]. This situation renders approaches based on scalar order parameters impractical. In what follows, we still name the low temperature state a gel phase, in spite of its ripple character.

The principle of the analysis is as follows. A 212 lipid molecules system was thermalized at low (288 K) and high temperature (358 K). At such temperatures, we assume that lipid conformations are predominantly gel and fluid respectively. Our training set was therefore composed of an equal number of conformations coming from the 288 K MD trajectory (all labelled as gel) and the 358 K trajectory (labelled as fluid). Details regarding the numerical simulation procedure and machine learning implementation details are given as Supplemental Material (SM).

Raw molecular conformations (spatial coordinates) were processed to remove all translation and rotation degeneracy, and the dimension of the initial conformation space was slightly reduced, resulting respectively in a 100-dimensional reduced coordinate space X , and a 61dimensional mutual distances space D. The coordinate space X and the distance space D provide a very detailed description of individual lipid conformations, and constitute the starting point of the ML approach. They are referred as "feature spaces" in the context of automated classification (see SI). The processed conformations were then fed to the ML algorithms. Three different algorithms were selected for the purpose of classifying the molecular conformations: Naive Bayes (NB), K-Nearest Neighbors (KNN) and Support Vector Machines (SVM). They were all used as implemented in the Python/Scikitlearn package [20]. The different ML algorithms were tested and found to perform moderately well when used separately, some methods performing better for gel lipids, other methods for fluid lipids (Fig S3 ,SM). Combining the approaches together, we managed to get a success rate of 88% upon validation, i.e. using 80% of the training set configurations for learning, our best ML algorithm was able to assign 86% of the low temperature configurations to a gel state, and 91% of the high temperature configurations to a fluid state. After training, the ML model was used to analyse simulations at arbitrary temperatures.

One important prediction of the two states model is the presence of minor components in the majority phase, under the form of "thermal excitations". A finite fractions of molecules tend to adopt a conformation different from their immediate environment, in spite of the presence of a "local field" biasing the statistics in favor of the dominant state. In that respect, the training set cannot be considered as containing only pure gel and fluid conformations. The presence of minor components is inherent to the presence of thermal excitations, and it does not seem possible to curate the training set without introducing further unwanted biases into the analysis. Our results show, however, that the training algorithm is not sensitive to the presence of a small fraction of non representative lipid conformations. In other words, the learning procedure was found to be robust so long as the training sets temperatures were chosen far apart from the melting transition.

Let us first analyze how ML predictions differ from those based on standard scalar structural order parameters. Figure (2) shows the distribution of the molecular segmental order parameter S mol , below and above the transition. This scalar observable was defined, for a given lipid and an instantaneous configuration (coming from a MD trajectory frame), by averaging over all the CC bonds in the two aliphatic tails a nematic parameter (3 cos(θ) 2 -1)/2, θ being the bond orientation with respect to the bilayer normal direction z. Both histograms overlap significantly. Using for instance a threshold value S th = 0.55, where the estimated probability density functions (pdf) coincide, it was found that altogether 17% of lipids (1 in 6) were set to be assigned to the opposite state, either gel at 358 K, or fluid at 288 K. The gel S mol distribution appears to be strongly skewed, as 28% (1 in 3) of the lipid molecules end up in the fluid state at 288 K. The relatively large fraction of lipids with S mol ≤ 0.55 found at 288 K is in part due to the ripple structure of the bilayer. This clearly show that using S mol for categorizing lipid conformations gives, at best, unreliable results. A similar analysis was conducted for two other scalar determinants: the lipid head-to-tail extension L, and the area per lipid A l computed from a Voronoi tessellation of the 2d xy projection of lipid center of mass positions (assigning to each lipid a polygonal cell of given area). Both were found to perform worse than the segment orientation S mol approach (Fig S5 ,SM).

To provide further evidence in favor of our approach, we compared the ML predictions to the membrane structure as a function of temperature. Figure (3) shows on the same graph the average area per lipid A l , the average volume per lipid V l and the ML prediction for the fraction of lipid assigned to the fluid state. The area per lipid is defined as the projected membrane area divided by the number of lipids per leaflet, neglecting out-of-plane membrane undulations. The volume per lipid results from the Voronoi analysis of the lipid center of masses (Fig S1 ,SM). The structural values evolve monotonically and reversibly with temperature, from 288 to 358 K. The three curves superimposes very well, which shows that structural data supports the finding of our classification tool. Two snapshots of the membrane upper leaflet are also provided, at low and intermediate temperature. The low temperature membrane shows a small number of isolated small clusters of fluid lipids. A one to one fluid/gel ratio is observed at 318 K (51 ± 3% of lipids in the fluid phase), a temperature close to the experimental melting temperature (314 K). A snapshot of the upper leaflet at 318 K shows two large distinct domains, separated by a FIG. 2. Distribution of the molecular order parameters S mol obtained from MD trajectories at 288 and 358 K respectively (using 10 6 lipid conformations). The dashed black line shows the threshold value that should be set to best determine the internal lipid states. The fraction of the population which would be incorrectly classified by using this threshold is represented by the shaded area in gray.

smooth boundary. The fluid state ratio evolves smoothly and reversibly from an estimated value of 14 ± 2% at 288 K to 95 ± 1% at 358 K. The very strong correlation between fluid ratio x f and area per lipid A l can be naturally interpreted in the framework of the two-states model by assigning constant values A f and A g to the fluid and gel states, resulting in

A l = x f A f + (1 -x f )A g .
A similar conclusion could be reached when considering the volume per lipid V l .

The ML predictions agree well with the global structural properties of the membrane and we now consider the local correlation properties. The two-states model combines the internal spontaneous dynamics of each lipid with local interactions promoting cooperativity between neighbors. The interactions, known as J coupling in the Ising model context, are responsible for the sharp structural and thermodynamic changes with temperature, the emergence of local correlations, clustering and domain formation, such as seen in Fig 3 . Such couplings create a local field whose effect is to bias the lipid state in favor of the dominant local phase, according to the majority rule. Internal reversible gel ↔ fluid state transitions occur spontaneously, according to a non conserved parameter, or Glauber dynamics [START_REF] Glauber | Timedependent statistics of the ising model[END_REF]. In order to get insight into the local correlations and flip rates, we performed a systematic statistical neighbor analysis, and estimated the conditional conversion states.

Voronoi tessellations provide an operational method for deciding unambiguously which lipid pairs are nearest neighbors. Following a protocol described in the SI, we consider two consecutive MD frames and divide the lipid population into four subsets: (a) lipids categorized as fluid, which stay in the fluid state, (b) lipids categorized as fluid, switching to the gel state, (c) lipids in the gel state, remaining in the gel state and (d) lipids in the gel state switching to a fluid state. At T = 318 K the fluid and gel phases compete evenly and the fluid and gel populations are roughly equal. We count for each molecule, the number n g of gel neighbors and n l of fluid neighbors. We find in these conditions that the most typical environment of a lipid in fluid state (a and b) is n g = 2 ± 1 and n f = 7 ± 1, while for a lipid in gel state (c and d) n g = 7 ± 1 and n f = 1 ± 1. However, the internal lipid state fluctuate spontaneously, and clearly the dynamics is strongly influenced by the local environment, as demonstrated in Figure (4). Indeed, for gel to fluid (case d) and fluid to gel (case b) transitions, we notice that all lipids that are just about to switch are more likely to have an equal number of gel and fluid neighbors (typically n g = 3 ± 1 and n f = 5 ± 1). We conclude that lipids subject to internal state transitions are mostly located at the border between domains. The results of our neighbor analysis clearly support the idea that internal state dynamics is under the control of some local field.

As a summary, we trained a Machine Learning algorithm to classify phospholipid molecular conformations obtained by atomistic molecular dynamics simulations. Lipids were sorted into two classes, gel and fluid, according to their similarity with a reference (training) set of conformations originating from low and high temperature trajectories respectively. The efficiency of the ML approach is superior to simple schemes based on scalar order parameters, and deals successfully with the ripple structure at low temperature. The measured fraction of fluid/gel conformations correlates very well with the observed structural changes as temperature evolves. A finite fraction of the minor phase is always present which can be associated to thermal excitations in the framework of a two-states Ising model interpretation. Lipids with similar state tend to cluster into large domains, while spontaneous internal state conversions are more likely to occur at the boundary between domains. The local distribution of neighbors supports the concept of local field and nearest neighbor coupling. We overall conclude that our ML approach provides convincing evidence in favor of the two-states phenomenological model.

We foresee numerous applications of the present approach. A first straightforward extension concerns lipid binary mixtures [START_REF] Ipsen | Modelling the phase equilibria in two-components membranes of phospholipids with different acyl-chain lengths[END_REF] and liquid ordered phases caused by the presence of cholesterol [START_REF] Ipsen | Phase equilibria in phosphatidylcholine-cholesterol system[END_REF], as well as membranes of more complex composition [START_REF] Ingòlfsson | Lipid organization of the plasma membrane[END_REF]. We also anticipate that our ML approach will be useful to study the influence of membrane solutes that are known to influence the thermodynamics of melting in model membranes. This includes hydrophobic pollutants, e.g. pyrene [START_REF] Fraňovà | Can pyrene probes be used to measure lateral pressure profiles of lipid membranes? perspective through atomistic simulations[END_REF], carbohydrates [18,[START_REF] Crowe | Solution effects on the thermotropic phase transition of unilamellar liposomes[END_REF] and synthetic oligomers [START_REF] Rossi | Polystyrene nanoparticles perturb lipid membranes[END_REF][START_REF] Rossi | Modeling the effect of nano-sized polymer particles on the properties of lipid membranes[END_REF]. A ML approach could then quantify the lipid state alteration induced by these compounds. Importantly, this analysis is also well-suited to study the local lipid environment of membrane proteins, for which the existence of lipid mediated interactions and minor phase nucleation is speculated [START_REF] Gil | Wetting and capillary condensation as means of protein organization in membranes[END_REF][START_REF] Gil | Theoretical analysis of protein organization in lipid membranes[END_REF]. For all these cases, significant improvements over approaches relying on scalar order parameters can be anticipated.

V.W. warmly thanks Tiago Espinosa de Oliveira for helping with simulations set-up. The authors gratefully acknowledge support from the high performance cluster (HPC) Equip@Meso from the University of Strasbourg, through grant n • G2018A53. directed along u z and the vector -----------→ C (j,k-1) C (j,k+1) defined by the carbon atoms surrounding atom k within the same tail j = 1, 2. The order parameter S mol,(j,k) of the atom k is obtained from the 2 nd Legendre polynomial P 2 using cos(θ (j,k) ) = u z • -----------→ C (j,k-1) C (j,k+1) , and averaging over j and k:

S mol = 1 2(N C -2) 2 j=1 N C -1 k=2 1 2 3 cos(θ (j,k) ) 2 -1 (1)

B. Next-nearest neighbors statistics

After completion of the 3d Voronoi tessellation using Voro++, a list of next-nearest neighbors was established for each lipid center of mass. We also collected the areas of the polygonal surfaces separating each pair of neighboring Voronoi cells. The ghost lipids and their corresponding faces were removed from the lists. The neighbor lists were further curated by removing all the neighbor pairs for which the corresponding face area accounted for less than 1% of the total surface area of each Voronoi cell in contact. The number of next-nearest neighbors were finally counted to build the coordination statistics (n g , n f ), where each lipid molecule has n g gel and n f fluid neighbors.

III. MACHINE LEARNING DISCRIMINATION OF LIPID STATES A. Lipid classification procedure

The lipid classification process involves three steps:

• The molecular conformation of a lipid is recorded as a list of atom positions. Lipids are then shifted and rotated in order to remove all rotation and translation degeneracy.

• The lipid molecular conformations are further simplified, reducing each single lipid conformation to a set of 100 (r, z) coordinates (configuration space X ), or 61 mutual distances (configuration space D, see below for details).

• 424 conformations were selected with the purpose of training the algorithms. A number of Machine Learning procedures were tested and compared. A combination of 4 algorithms was found to maximize the training success rate. We therefore combine the 4 algorithms in our final classification procedure (see below for details).

B. Definition of the configuration spaces

Machine Learning algorithms share the ability to discriminate well multidimensional data. Our purpose is to feed the algorithms with single lipid conformations and have them sorted into two classes. We observe that the raw chemical formula of DPPC (CAS 2644-64-6) is PNO 8 C 40 H 80 . Disregarding hydrogens, each lipid molecule comprises 50 "heavy atoms", and therefore the associated single lipid configurations belong to a 150 dimensions vector space. As some configurations can be mapped onto each other by means of a spatial displacement (translation and/or rotation), the set of configurations has only 144 independent degrees of freedom.

For simplicity, we decided to work with slightly smaller configuration spaces, by projecting further the 144 dimensional molecular conformations onto two configuration spaces, X and D, that we now define. Within the first approach, one determines the proper inertial frame of each lipid configuration, locating the center of mass, computing and diagonalizing the matrix of inertia. Lipid configurations can be recast into the inertial frame, using the center of mass as origin, and the longest axis (smaller moment of inertia) as vertical axis. This almost always results in having the lipid directed along the bilayer normal, with two possible orientations. When necessary, the lipid is flipped in order to place by convention the phosphocholine group into the positive z upper half space. This combination of reorientations lifts entirely the translation and orientation degeneracy of the original coordinate space. New cylindrical coordinates {r i , θ i , z i }, i = 1 . . . 50 can be associated to the resulting lipid conformation. We decided then to disregard entirely the angles θ i , keeping only the coordinate subset {X i = r i , z i }, i = 1 . . . 50. This defines a N X = 100 dimensional space X that will subsequently be referred as "Coordinate space".

For the second approach, starting from the original 150 lipid coordinates, we calculate a set of mutual euclidean distances between pairs of atoms. Two atoms participate in a pair when they are separated by 6 positions along In both cases, the original dimension of the problem is reduced (from 144 to 100 and 61 respectively) but yet the configuration spaces X and D are large and preserve to a large extent the complexity of the original conformations. As our comparisons show, the classification of lipid states is efficient, whether one starts from X or D. Three different algorithms commonly used for Machine Learning classifications were used in this Letter. A short description of these algorithms is given below.

C. Naive Bayes

Given a configuration space D = {D j }, j = 1 . . . N D , and two classes s = {gel, fluid}, the Bayesian approach assumes that it exists a joint probability distribution P (s, {D j }), that weighs the respective likelihood of each state s once a configuration {D j } is provided. The Bayesian model makes a decision regarding the state s by comparing the conditional probability densities P (gel|{D j }) and P (fluid|{D j }). To be precise, the Bayesian approach attributes a fluid label to a configuration if the sign of ln P (fluid|{D j })

P (gel|{D j }) (2)
is positive, and a gel label otherwise.

Each Bayesian approach provides a mathematical model P ({D j }|s) describing the expected configuration distribution for a given class s: fluid and gel. There is in principle entire freedom in choosing the model, but the efficiency and the optimization requirements limit such choices in practice. Training a Bayesian model therefore amounts to finding the most realistic function P ({D j }|s) as far as classifying a given training set of data is concerned.

The Bayes theorem provides the connection between the conditional probabilities entering in the choice function (2) and the model:

P (s|{D j }) = P ({D j }|s)P (s) P ({D j }) , (3)
where P (s) represents the prior distribution of s, i.e. the statistical distribution of s in the absence of any configuration related information, and P ({D j }) = P (gel, {D j }) + P (fluid, {D j }) a normalization factor which cancels out in eq (2). The Naive Bayes (NB) gaussian model assumes that P ({D j }|s) factorizes as a product of independent gaussian distributions of D j .

P ({D j }|s)

N D j=1 dD j = N D j=1 dD j √ 2πσ s,j exp - (D j -D s,j ) 2 2σ 2 s,j (4)
Training the algorithm means finding the best mean value D s,j and standard deviation σ s,j for every parameter D j in the distance space D and each class s. The number of parameters to compute turns out to be equal to twice the dimension of D. In our case, the training set contains an equal number of gel and fluid conformations, and there is no a priori bias between classes, meaning that P (gel) = P (fluid) = 1/2. Therefore, the number of parameters to determine during training is 2 × 61 = 122. To sum up, the Naive Bayes approach classifies a lipid conformation by computing a quadratic score function in the conformation space,

N D j=1 (D i -D fluid,i ) 2 2σ 2 fluid,i - N D j=1 (D i -D gel,i ) 2σ 2 gel,i + Const (5)
and deciding whether a {D i } lies closer to {D fluid,i } or to {D gel,i } according to this generalized distance.

D. K-Nearest Neighbors

The K-Nearest Neighbors (KNN) classification method is based on defining a distance between any arbitrary pairs of objects to discriminate. A natural choice is the Euclidean norm of the feature space, here the N X dimensional coordinate space X .

The KNN algorithm finds the K nearest neighbors within the training set, of each new configuration to classify. Decision is taken based on majority rule, i.e. the most abundant class found among the K closest neighbors. The optimal K is determined during the training and validation process, and was set equal to 5, a typical value for this method.

E. Support Vector Machines

Support Vector Machines (SVM) classify data by means of linear separation in high dimensional representation spaces. Denoting φ an arbitrary element in a given representation space R, the binary decision is given by the sign of the affine expression w • φ + b, with b a numerical constant and w the hyperplane of separation normal vector. In a few favorable cases, it is possible to use directly the data definition space as representation space. However, in many practical situations, efficient classification can only be achieved by mapping the data (e.g. X or D) onto a larger representation space, namely x → φ(x). Training a SVM corresponds to choosing a suitable representation space R, and finding the optimal b and w. As shown in [18,19], given a training set {x i }, i = 1 . . . n, the optimal w can always be expressed as a linear combination w = i α i y i φ(x i ) with either positive or vanishing α i coefficients, and y i = ±1, depending on the class (fluid 1, gel -1) of the corresponding vector data x i .

The subset of vectors {x i } participating in the definition of w with non vanishing coefficients α i > 0 forms the so-called support vectors. Denoting J the sequence of indices of support vectors, and K(x, x ) the product φ(x)•φ(x ) in R, called kernel function, the SVM decision function for any vector data x reads:

sign   b + j∈J α j y j K(x j , x)   . (6)
The SVM training optimization problem can therefore be formulated without any explicit reference to the representation space R, nor the mapping φ(x). It only requires an explicit positive kernel function K(x, x ). As explained in [18,19], there are efficient quadratic optimization algorithms for determining the support vectors x j , the non-vanishing coefficients α j and the shift constant b.

In this study, we used the standard radial basis kernel function

K(x, x ) = exp -γ||x -x || 2 , (7)
with a default value γ equal to the inverse of the dimension of the configuration space. This choice assumes that each component of x is of order 1. When using the coordinate space X , the trained SVM ends up using 297 non vanishing support vectors and coefficients α j , γ = 1/100, b = -0.416 (167 fluid, 130 gel support vectors). When considering the distance space D, the trained SVM used 161 non vanishing coefficients α j , γ = 1/61 and b = 0.134 (100 fluid, 61 gel support vectors).

IV. TRAINING

The Machine Learning analysis performed in this Letter were conducted using the Scikit-Learn module for Python 3 [20]. An unbiased selection of lipid conformations extracted from a trajectory at low temperature (288 K) was part

V. ASSERTING THE PREDICTIVE CAPACITY OF EACH MODEL

Once defined the training set, with properly labelled gel and fluid states, 20% of the lipids in each phase were taken apart for forming a validation set. After training the models (i.e. optimizing the parameters with respect to the 80% remaining conformations) a prediction score was separately calculated for the gel and fluid conformations in the validation set. The overall procedure was repeated 10 times, each time with the same training set, but independently drawn validation subsets.

The resulting scores shown in Fig. S3 shows that the NB predictive capacity is independent from the phase of the lipid considered. KNN performs better for fluids than for gels. At the opposite, "coordinates" X -SVM and "distances" D-SVM seems to perform better in the gel phase.

We therefore decided to combine the predictions of the above models, retaining those who perform best in each phase. The following decision chain was implemented:

1. If the 4 models agree on the same prediction for the phase, this prediction is retained; 2. if the D-SVM algorithm predicts a gel phase, the lipid configuration is assumed to be gel ; 3. if the X -SVM algorithm predicts a fluid phase, the lipid configuration is assumed to be fluid ; 4. if the NB algorithm predicts a gel phase, the lipid configuration is assumed to be gel ; 5. if none of the above conditions have been met, the KNN algorithm makes the final decision on the configuration classification. The order parameter curves (Fig. S4(a)) clearly discriminate among the low temperature (288 K) and high temperature (358 K) fluid phases, in agreement with published results on these systems [21].

The order parameter of the atoms in the lipid tails at low (288 K) and high (358 K) temperature is characteristic from membranes in the gel and fluid phases respectively [22]. The phase transition can be clearly seen in Fig. S4(b) as a significant variation in the evolution of the area per lipid A l around 321 K.

Experimental structural values are available at 323 K [23][24][25][26]. Nagle et al. obtained for DPPC an area per lipid equals to 64 ± 1 Å2 significantly close to the value circa 60 Å2 we obtained in our simulations. Using the average DPPC bilayer thickness reported by Nagle et al., we could estimate an experimental volume per lipid of 1220 ± 50 Å3 , which agrees fairly with our Voronoi value of 1300 Å3

VII. NAIVE CLASSIFICATIONS

The distributions of the areas per lipid and molecular elongations at low and high temperatures are shown in Fig. S5. Using a naive classification scheme based on a single threshold value for either of the two previous scalar parameters would at best result in a prediction accuracy of respectively 69% and 67%.

FIG. 1 .

 1 FIG. 1. Snapshots of a DPPC bilayer simulated at (top-left) 288 K and at (bottom-left) 358 K, respectively below and above the experimental Tm of the lipid. Lipid molecules are shown in gray, with their phosphorus atom displayed as a plain big sphere to distinguish the lipid orientation in the bilayer. Blue lines delimit the simulation box beyond which periodic boundary conditions are applied to the system. Right: lipids extracted from the bilayers displaying an important diversity in their conformations, with their associated molecular order parameter S mol . Top-right: conformations at 288K, bottom right: conformations at 358K.

FIG. 3 .

 3 FIG. 3. Ratio of lipids in fluid state (black curve, vertical axis on the left), area per lipid (blue curve, vertical axis on the right) and volume per lipid (red curve, vertical axis on the right) as function of temperature. Two snapshots of the upper membrane leaflet, respectively taken at 288 and 318 K, are shown in the inset. The dark spots correspond to the 2d projection of the lipid center of masses. Boundaries between lipids result from the associated 2d Voronoi tessellation. Each cell color corresponds to the assignment of the ML, blue for the gel state and red for the fluid state.

FIG. 4 .

 4 FIG. 4. Color coded histograms of the (n l , ng) distribution for the 4 subsets described in the text, at T = 318 K. Probability at T = 318K for a lipid to have its local environment composed of ng neighbors in gel state and n f neighbors in fluid state. Left column informs about the local environment probability of lipids in the fluid state (Top) and gel state (Bottom) when no change in the state is observed. Right column shows the probability of the local environment just before lipids experience a transition to the other state.

  FIG. S1. Comparison of 3-dimensional Voronoi tessellations of a lipid bilayer configuration (A) without and (B) with ghost lipids. Without ghost lipids, most cells are unbounded, with infinite volume, due to the absence of particle on the opposite side of the water-membrane interface. As a practical solution of this problem, ghost lipids are added to the data set, mirroring the lipid center of mass positions. The resulting cells for lipids inside the bilayer display a realistic volume and shape, accounting for the water interface in a natural way.

  FIG. S2. Skeletal structure of a DPPC molecule. The atoms are colored according to the number of pairs formed between this atom and other atoms separated by exactly 6 covalent bonds along the molecule chain. The highest number of neighbors for a single atom was found to be 6. The total number of pairs is 61.

  FIG. S3. Prediction scores of the Machine Learning classification methods for lipids in the gel phase (Left) and in the fluid phase (Right). Besides the Naive Bayes method (NB), all methods have an important asymmetry in their accuracy between each phase. We call these asymmetries the "expertises" of the methods.

  FIG.S4. Confirmation of the thermodynamic phase of the bilayer using two common structural parameters: the order parameter S mol and the area per lipid A l . (Left) the average order parameter is shown as a function of the carbon atom index along the chain (from glycerol to terminal end), for each sn1 and sn2 chain. (Right) the phase transition can be clearly seen in the evolution of the area per lipid as a function of temperature. A sigmoid fit points to a transition temperature Tm equal to 321 K in our system.

Fig

  Fig.S6compares the histogram of molecular order parameters S mol as a function of the temperature of the lipid bilayer from which the configurations are extracted (288 K or 358 K), and as a function of the result of the Machine Learning classification procedure. The difference between the distribution at 288 K and the distribution in the a posteriori gel state ensemble indicates that a small fraction of lipids in the fluid state are already present at 288 K.

  FIG. S5. Distributions of the area per lipid (Left) and the average elongation between phosphorus and sn1 terminal carbon atoms (Right) from lipids conformations at 288 and 358 K.

Supplementary Information

I. SIMULATIONS

All simulations were performed using GROMACS 2016.4 [1,2] along with the CHARMM-36 all-atom force-field [3] (June 2015 version). A lipid bilayer made of 106 lipid molecules per leaflet, each containing 130 explicit atoms, was created using CHARMM-GUI [4][5][6][7]. It was hydrated with two 8 nm thick water layers on each side (connected through periodic boundary conditions), using the TIP3P water model, for a total of 29826 solvent molecules. The force field parameters for DPPC molecules were provided directly by CHARMM-GUI [8,9].

The above system was first subject to energy relaxation using steepest descent energy minimization, followed by a 10 ps NVT thermalization stage at 288 K. Then, the bilayer was subject to a 1 ns NPT equilibration run coupled to a semi-isotropic barostat (1 bar in all directions). The system was then further equilibrated at the desired temperature (pressure conditions unchanged) during a second NPT equilibration step of 10 ns. Molecular dynamics production runs of 50 ns were finally generated at the same temperature and pressure conditions. The analysis were performed on the last 25 ns of simulations. All time steps were set to 2 fs.

All the molecular dynamics simulations used the leap-frog integration algorithm [10]. Temperature and pressure were kept constant using respectively a Nosé-Hoover thermostat [11,12] (correlation time τ T = 0.4 ps) and a Parrinello-Rahman barostat [13,14] (correlation time τ P = 2.0 ps, compressibility 4.5 × 10 -5 bar -1 ).

Lipid and water molecules were separately coupled to the thermostat. Following GROMACS recommendations for the CHARMM-36 all-atom force field, a Verlet cut-off scheme on grid cells was used with a distance of 1.2 nm, and non-bonded interactions cut-offs (Van der Waals and Coulombic) were also set to 1.2 nm. Fast smooth Particle-Mesh Ewald electrostatics was selected for handling the Coulombic interactions, with a grid spacing of 4 nm. A standard cut-off scheme with a force-switch smooth modifier at 1.0 nm was applied to the Van der Waals interactions. We did not account for long range energy and pressure corrections, and constrained all the hydrogen bonds of the system using the LINCS algorithm.

II. SYSTEM ANALYSIS A. Determination of structural parameters

Values of the average area per lipid A l and order parameter S mol of the bilayer were respectively obtained using the GROMACS built-in commands gmx energy and gmx order.

For measurements of individual lipid properties (area, order parameter, elongation), atom positions were collected from trajectories using the Python 3 MDAnalysis module [15,16]. The individual areas per lipid were obtained from Voronoi tessellations of the two-dimensional projections of the lipid center of masses, computed using the Voro++ library [17]. The individual volumes per lipid were derived from three-dimensional tessellations of the lipid centers of mass, again using the Voro++ library. Note that the bilayer geometry requires a specific tessellation procedure: this was done by introducing ghosts lipids in the water regions. Without these ghost lipids, the tessellation cells cannot be correctly defined and are unbounded across the membrane-water interface, thus overestimating significantly the individual volume per lipid. Ghosts lipids are mirror images of bilayer lipids across the local lipid-water interface (cf. Fig. S1). After the tessellation was made, ghost lipids and their corresponding cells were removed, and only the volumes of physical lipids were collected and analyzed.

The molecular order parameter S mol of individual lipids was calculated by measuring, for every N C -2 = 14 non-terminal carbon atoms k = 2 . . . 15 of the 2 tails of the lipids, the angle formed between the z-axis of the system


Tiago Espinosa De Oliveira

Fabien Leonforte

Nicolas-Morgantini Luc

Anne-Laure Fameau

Bernard Querleux

Fabrice Thalmann

Carlos M Marques

A fluid bilayer phase in aqueous mixtures of fatty alcohol and cationic surfactant

Lamellar gel networks are creamy formulations extensively used in cosmetics and pharmaceutics.

The basic units of these networks are lipid bilayers assembled from mixtures of fatty alcohols, surfactants and water. Despite being an essential step for the preparation of the dispersions, the nature of the mixtures at high temperature has defied understanding. Here, we show by atomistic molecular dynamics simulations that aqueous mixtures of cetyl (C 16 OH) and stearyl (C 18 OH) alcohols, and cetyl-trimethylammonium chloride (CTAC) lead to fluid bilayers, the units of the lamellar L α phase. By providing a consistent numerical simulation model able to describe the structure and properties of fatty alcohol lamellae stabilized by surfactant, our work paves the way for the elucidation of the forces regulating this family of industrially important gels.

Fatty alcohols (FA) are amphiphilic molecules associating a polar hydroxyl head to an alkyl chain tail (CH 3 -(CH 2 ) n-1 -OH) [START_REF] Noweck | Ullmann's Encyclopedia of Industrial Chemistry[END_REF]. Cetyl (n = 16) and stearyl (n = 18) alcohols, in combination with non-ionic or ionic surfactants, are common ingredients in formulations of cosmetic and pharmaceutical cream products [START_REF] Barry | [END_REF][3][4][5][6][7]. These aqueous dispersions, often known as lamellar gel networks, derive their advantageous properties from an extended and highly interconnected lamellar structure [START_REF] Iwata | Cosmetic science and technology: theoretical principles and applications[END_REF]. The basic network units are rigid (L β ) bilayers, or stacks of bilayers, forming a percolating structure that can withstand elastic deformations [START_REF] Benton | [END_REF][10][11].

Depending on temperature, fatty alcohols alone adopt various ordered structures with poor hydration capacity [12][13][14], preventing mixtures of water and fatty alcohols to be formulated as creams with the desired homogeneity, stability and viscoelasticity: only when an appropriate proportion of surfactants is added to the mixture can the lamellar gel network structure be obtained [START_REF] Iwata | Cosmetic science and technology: theoretical principles and applications[END_REF].

The textural properties of the cream are empirically optimized by various preparation steps involving combinations of stirring and heating or cooling. But in all formulation variants the different components are first brought together at high enough temperatures where a fluid bilayer (L α ) phase is believed to form [15]. Despite being a key determinant for the formation of the lamellar network, the nature of such high temperature phase of fatty alcohols/surfactant mixtures has not been throughly inspected. In particular, to our knowledge there are no published molecular dynamics simulations of mixtures of water, fatty alcohol and surfactants, thus precluding the long-sought elucidation of the forces shaping these membranes and controlling their interactions.

In this Letter we perform all-atom molecular dynamics simulations of aqueous solutions of cetyl (C 16 OH) and stearyl (C 18 OH) alcohols and cetyl-trimethylammonium chloride (CTAC), a typical mixture leading to the formation of lamellar gel networks [START_REF] Iwata | Cosmetic science and technology: theoretical principles and applications[END_REF]10,15].

Cetyltrimethyl-ammonium chloride (CTAC) is a cationic surfactant that mixes well with cetyl and stearyl alcohols that have similar chain lengths. Fatty alcohols (FA) and surfactant (CTAC) molecules were parameterized using the CHARMM 36 [16,17] force-field.

In particular, we used for the alkyl chains the same parameters as the ones introduced initially for modeling 1,2-dipalmitoyl-phosphatidylcholine (DPPC) tails. As a matter of fact, CHARMM-36 was shown to accurately describe the melting transition of DPPC lipid bilayers, as a result of some empirical optimization. Methanol CH 3 OH, and DPPC phospho-choline parameters were used respectively to create the alcohol and trimethyl ammonium groups. The TIP3P water model was used for the solvent.

Most of the simulated systems were made of 256 chains arranged into a symmetric bilayer conformation, each opposing leaflet containing 128 chains, and solvated with 2560 water molecules. Initial configurations were prepared by molecular modeling. Four chemical compositions were considered: pure cetyl (C 16 OH), 20% molar surfactant (C 16 OH:CTAC (4:1)), equimolar cetyl-stearyl mixtures without (C 16:18 OH) and with 20% molar (C 16:18 OH:CTAC (4:1)) surfactant. Triplicate 100 ns runs for annealing and quenching were performed to thermalize the systems [18]. This choice for initial configurations was justified by a number of available experimental results consistent with a lowly hydrated gel L β structure at room temperature [15].

Structural changes in the FA organisation were followed by computing the carbon-carbon order parameter S of the alkyl chain with respect to the bilayer normal direction z. It is defined for every group (CH 2 ) i in a chain as:

S i = 1 2 3cos 2 θ -1 , (1)
where θ is the angle between the two nearest carbons surrounding the selected CH 2 group and the membrane normal axis. Bilayer systems are characterized by a positive S, which is larger for gel than for fluid phases. Isotropic structures lead to a vanishing order parameter S i = 0. In all the cases of interest, the order parameter decreases as i shifts from the FA head towards the terminal methyl group (-CH 3 ). An average order parameter S, representative from the central section of FA chains (averaged over 8 carbons, numbered from 4 to 11), was computed and represented in Fig. 1.

In the absence of surfactant, both C 16 OH and C 16:18 OH showed a well defined gel phase at temperatures below respectively 298 and 308 K, with S ∼ 0.7. In this phase, alkyl chains adopt well packed dense all-trans conformations. At larger temperatures, these two systems were found in a completely disordered phase, presumably an isotropic fluid phase. The preferential arrangement of FA molecules into leaflets was lost. Cetyl-stearyl mixtures melted at a slightly higher temperatures than pure cetyl systems, in agreement with experimental data [19].

In the presence of a 1:4 CTAC surfactant ratio (20% of the system molar mass), a structure consistent with a lamellar fluid phase was observed in a [323-343 K] temperature range. Tails were more disordered, melted, as shown by the lower value S ∼ 0.35. Experimental systems of comparable FA-CTAC composition but with much larger water contents display a L β → L α transition temperature at 243 K. In our simulations, above 358 K, an isotropic fluid phase, similar to the one seen without CTAC prevailed. We observed that repeated cycles of annealing-quenching resulted in a denser, better packed low temperature solid structure, with an average order parameter higher than 0.7 -see Fig. 2.

As further evidence of a lamellar fluid state stabilized by surfactants, we increased the size of the simulated systems and tested the ability of the bilayer to spontaneously form a vesicle, a topological hallmark of fluid bilayers. A large system comprising 2560 fatty chains was initially prepared as a flat bilayer. Then, a flat square slab was cut off and dunked into a large water reservoir, in such a way that the slab did not extend across periodic boundary conditions (PBC) and had a free boundary line exposed to the solvent (Fig. 3 at 0 ns). As the time sequence of snapshots displayed in Fig. 3 shows, the slab spontaneously closed on itself, adopting a small unilamellar vesicle shape. From 5 ns to 20 ns, the bilayer formed a temporary bicelle disk, closed by a folded leaflet rim. This disk-like structure remained for some time, while undergoing severe surface fluctuations. Then, it started to bend spontaneously after 30 ns and ended up forming the closed vesicle after 40 ns of simulation time. The vesicle stayed stable during the following 500 ns of simulation.

The lamellar fluid phase structure was characterized further by measuring the area per fatty chain A f , and the membrane thickness D m . Both parameters were found to evolve smoothly with temperature except for an abrupt change between 315 and 323 K, supporting evidence of a sharp gel-fluid transition (Table I). The bilayer isothermal stretching elastic modulus K A was estimated based on equilibrium box area fluctuations [20][21][22][23][24], according to the relation:

K A = k B T A A 2 -A 2 , ( 2
)
where A and A 2 are the average value and the mean square fluctuation (variance) of the sample monolayer area A = L x L y , k B the Boltzmann constant and T the temperature.

The area fluctuations arise from coupling the system to a semi-isotropic Parrinello-Rahman barostat [25] and a Nose-Hoover thermostat [26,27]. The Parrinello-Rahman barostat is assumed to enforce a constant pressure and vanishing surface tension NPγT ensemble, i.e. by allowing the L z and L x = L y lateral box sizes to be independently rescaled. We found a significant drop in K A at the transition (Table I), consistent with a lamellar fluid L α phase, with K A values comparable to those of lipid bilayers (∼ 250 mN/m) [28]. Note that imposing a semi-isotropic barostat is only possible if the simulated system possesses a finite stretching modulus K A . The highest temperature states occurring in the presence and in the absence of CTAC were shown to be unstable under a semi-isotropic barostat, displaying unbound L x -L z fluctuations, as expected from isotropic fluid phases. Gel-fluid bilayer transitions are commonly believed to be weakly first order [START_REF] Heimburg | Thermal Biophysics of Membranes[END_REF]. To gain insight into the thermodynamic transition properties, we estimated the change in enthalpy ∆H upon melting. The enthalpy of the bilayer and water system H was defined as H = U + P V [START_REF] Seo | [END_REF], where U is the time averaged total energy, computed from the sum of all the potential energies in the force field and the translation kinetic energy of the atoms, temperature phases. We found respectively for these jumps a value of 10.30 and 10.00 kJ/mol. This result is quantitatively consistent with experimental differential scanning calorimetry data (DSC) [15], and falls within the expected range of chain melting values for equivalent phospholipid systems [START_REF] Marsh | Handbook of Lipid Bilayers[END_REF]. We found that the surfactant containing cetyl and cetyl-stearyl mixtures melts at different temperatures, the longer chain alcohol driving L β → L α melting to higher temperatures, albeit with comparable transition enthalpies.

T (K) A f (nm 2 ) D m (nm) K A (mN/
P
To summarize, we have proposed an atomistic model that successfully describes the cetyl and cetyl-stearyl alcohols mixtures in aqueous solution with and without surfactants (CTAC). We investigated the structural, mechanical and thermodynamic properties of these systems at various temperatures. In the absence of surfactant, we found a transition between a solid gel bilayer and an isotropic fluid phase in good agreement with experiments.

In the presence of surfactant (20% CTAC), we found 2 transitions. A first transition was seen between a gel phase and fluid lamellar phase. Our estimate of the enthalpy of melting

FIG. 1 .

 1 FIG. 1. Average order parameter as a function of temperature for four different compositions. (a) Aqueous C 16 OH and C 16 OH:CTAC (4:1) and (b) aqueous C 16:18 OH and C 16:18 OH:CTAC (4:1).

FIG. 3 .

 3 FIG. 3. Spontaneous formation of a small unilamellar vesicle from a C 16:18 OH:CTAC (4:1) fatty alcool-surfactant mixture at 333 K.

  (4:1): A f mean crosssectional area per chain, D m bilayer thickness. The relative statistical uncertainty for A f and D m is about 0.003. Bilayer elastic coefficient K A , from a fluctuation analysis argument. The relative statistical uncertainty turns out to be larger, of the order of 20% to 30%. All three parameters change sharply between 315 and 323 K. Parameters of the bilayer are measured over the last 80 ns of 100 ns simulations.

FIG. 4 .

 4 FIG. 4.Enthalpy as a function of temperature for the C 16 OH:CTAC (4:1) and C 16:18 OH:CTAC (4:1) systems.

Figure 4

 4 Figure 4 represents the enthalpy difference H(T ) -H(278 K) as a function of temperature for C 16 OH:CTAC (4:1) and C 16:18 OH:CTAC (4:1) aqueous solutions. The enthalpy curve shows a smooth trend (constant pressure specific heat) except in the region of the melting transition. The enthalpy variation at the transition of C 16 OH:CTAC (4:1) and C 16:18 OH:CTAC (4:1) was estimated by extrapolating the linear trends of the low and high

TABLE I .

 I Geometrical parameters for gel and fluid bilayers of C 16 OH:CTAC

				m)
	278	0.23	4.07	2328.4
	288	0.23	4.08	2274.1
	298	0.23	4.09	2146.3
	308	0.23	4.10	1875.4
	315	0.23	4.20	2077.2
	323	0.29	3.49	216.9
	333	0.30	3.47	226.9
	343	0.31	3.43	173.9

is very satisfactory although, perhaps not surprisingly, our temperature transition is lower than that experimentally measured in more hydrated multilamellar systems [15]. We found also that the compressibility modulus falls abruptly upon melting to a value comparable to similar known L α phases. The parameters used in the present work describe well the systems under investigation, both from the thermodynamic and from the mechanic point of view, and we confidently associate the lamellar fluid phase to a L α phase, while the low temperature gel phase is consistent with a L β state. Simulation of larger system sizes showed that the fluid bilayer may form stable unilamellar vesicles. This work is a first step into a better assessment of the influence of the high temperature phase of mixtures of fatty alcohols and surfactants on the formation of lamellar gel networks.

We anticipate that it will allow for instance exploring the role of the surfactant molar fraction in the stabilization of the fluid bilayer phase and to draw phase diagrams that can be directly compared to experiments. Importantly, the force-field model presented here paves the road for extending MD simulations into multi-bilayer systems, where one can test the influence of the bilayer-bilayer interactions on the build-up of the lamellar network . T.E.O., F.L., L.N.M., A.L.F., B.Q., F.T., and C.M. acknowledge partial support from the Investissements d'Avenir program "Développement de l' Économie Numérique" through the SMICE project. T.E.O., F.T., and C.M. acknowledge partial support from L'Oréal Research & Innovation, and F.L., L.N.M., A.L.F., and B.Q. are employees of L'Oréal Research & Innovation. We acknowledge support from the HPC cluster TGCC (Très Grand Centre de Calcul du CEA, Bruyères-le-Châtel, France).


Kaarel Mänd
email: kaarel.mand@ualberta.ca.

Stefan V Lalonde

Leslie J Robbins

Marie Thoby

Kärt Paiste

Timmu Kreitsmann

Päärn Paiste

Alexandr Christopher T Reinhard

Alexander E Romashkin

Noah J Planavsky

Kalle Kirsimäe

Aivo Lepland

Kurt O Konhauser

Palaeoproterozoic oxygenated oceans following the Lomagundi-Jatuli Event

come

The ~2220-2060 Ma Lomagundi-Jatuli Event (LJE) was the longest positive carbon isotope excursion in Earth history and is traditionally interpreted to reflect increased organic carbon burial and a transient rise in atmospheric O2. However, it is widely held that O2 levels collapsed for more than a billion years following the LJE. Here, we show that black shales postdating the LJE from the ~2000 Ma Zaonega Formation (ZF) contain the highest redox-sensitive trace metal concentrations reported in sediments deposited before the Neoproterozoic. This unit also contains the most positive shale U isotope (δ 238 U) values measured to date, providing novel evidence that there was a transition to modern-like biogeochemical cycling during the Paleoproterozoic. Although these records do not preclude a return to anoxia during the Paleoproterozoic, they uniquely suggest that the oceans remained well-oxygenated millions of years after the termination of the LJE.

The Paleoproterozoic era (2500-1600 Ma) witnessed the longest positive carbon (δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb) isotope excursion in Earth's history, recorded in marine carbonates deposited worldwide between ~2220-2060 Ma-the Lomagundi-Jatuli Event (LJE) [START_REF] Schidlowski | Carbon isotope geochemistry of the Precambrian Lomagundi carbonate province, Rhodesia[END_REF]2 . Marine carbonates are generally characterized by δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb values that vary between -5‰ to 5‰ for most of Earth's history, whereas peak LJE δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb values reach 10-15‰ [START_REF] Schidlowski | Carbon isotope geochemistry of the Precambrian Lomagundi carbonate province, Rhodesia[END_REF][START_REF] Karhu | Carbon isotopes and the rise of atmospheric oxygen[END_REF][START_REF] Bekker | Oxygen overshoot and recovery during the early Paleoproterozoic[END_REF] . There is ongoing debate regarding the processes that led to the LJE. The standard and most commonly accepted interpretation invokes increased burial of [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] C-depleted organic carbon (Corg), leading to the [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] C-enrichment of the dissolved inorganic carbon pool and positive δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb 3,[START_REF] Bekker | Oxygen overshoot and recovery during the early Paleoproterozoic[END_REF] . This acceleration in Corg burial may have been the result of an increase in oxidative weathering and nutrient delivery to the oceans [START_REF] Bekker | Oxygen overshoot and recovery during the early Paleoproterozoic[END_REF] . Using a simple isotope mass balance, it has been estimated that between 5-9×10 20 mol of Corg was buried over 100 Ma, corresponding to the release of O2 equivalents representing 12-22 times the present O2 atmospheric pool. An ensemble of evidence points to atmospheric O2 accumulation during the LJE, including proxies indicating growth of the marine sulphate reservoir [START_REF] Planavsky | Sulfur record of rising and falling marine oxygen and sulfate levels during the Lomagundi event[END_REF][START_REF] Scott | Pyrite multiple-sulfur isotope evidence for rapid expansion and contraction of the early Paleoproterozoic seawater sulfate reservoir[END_REF][START_REF] Blättler | Two-billion-year-old evaporites capture Earth's great oxidation[END_REF] , elevated concentrations of redox sensitive elements (RSE) [START_REF] Konhauser | Aerobic bacterial pyrite oxidation and acid rock drainage during the Great Oxidation Event[END_REF][START_REF] Partin | Large-scale fluctuations in Precambrian atmospheric and oceanic oxygen levels from the record of U in shales[END_REF][START_REF] Kipp | Selenium isotopes record extensive marine suboxia during the Great Oxidation Event[END_REF][START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] , and evidence for locally oxic conditions [START_REF] Bekker | Oxygen overshoot and recovery during the early Paleoproterozoic[END_REF][START_REF] Canfield | Oxygen dynamics in the aftermath of the Great Oxidation of Earth's atmosphere[END_REF] .

However, others have argued that the 'standard' interpretation of the LJE is difficult to reconcile with our understanding of the C and O cycles [START_REF] Bachan | The rise of oxygen and siderite oxidation during the Lomagundi Event[END_REF][START_REF] Shields | Tectonic controls on the long-term carbon isotope mass balance[END_REF][START_REF] Miyazaki | Making sense of massive carbon isotope excursions with an inverse carbon cycle model[END_REF] , and that there is a notable paucity of Corg-rich deposits at this time [START_REF] Kump | Isotopic evidence for massive oxidation of organic matter following the Great Oxidation Event[END_REF] , the most basic tracer of enhanced organic carbon burial. Further, the carbon isotope dynamics during and after the positive excursion are currently debated, and there are multiple interpretations of the LJE positive δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb values that do not invoke elevated rates of organic carbon burial and oxygen release (ref 18). Therefore, additional constraints on carbon and oxygen cycling are clearly needed to refine our view of this time interval.

In this work we present new redox tracer data from three correlative cores (OnZaP-1, 3 and OPH) of the upper Zaonega Formation (ZF), a mudstone-dolostone succession deposited in the aftermath of the LJE. We provide some of the most straightforward evidence for Proterozoic surface oxygenation. In fact, current redox proxy records are most consistent with increasing oxygenation in the aftermath of the LJ positive carbon isotope excursion-forcing a re-evaluation of our basic view of this turbulent interval of Earth's history.

Age and geochemical signatures of the Zaonega Formation

According to the latest U-Pb age constraints on tuff horizons within the lower ZF, the formation was deposited at <1970 Ma, ~70 Myrs after the proposed termination of the LJE [START_REF] Martin | Multiple Palaeoproterozoic carbon burial episodes and excursions[END_REF] . This age model is consistent with δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb stratigraphy. In the Onega Basin, δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb values in the lowermost ZF and underlying Tulomozero Formation are characterized by typical LJE values of ≥8‰ [START_REF] Melezhik | Carbonate deposition in the Palaeoproterozoic Onega basin from Fennoscandia: a spotlight on the transition from the Lomagundi-Jatuli to Shunga events[END_REF] , with δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb in preserved carbonate strata returning to normal marine values of ~0‰ further upsection [START_REF] Kreitsmann | Hydrothermal dedolomitisation of carbonate rocks of the Paleoproterozoic Zaonega Formation, NW Russia -Implications for the preservation of primary C isotope signals[END_REF] . In the OPH core, the studied section occurs after several hundred meters of stratigraphy-including cumulatively ~210 meters of mudstone-dominated sediments-and contains carbonates bearing a negative or normal marine (and therefore post-LJE) δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] Ccarb signal (ref. 30). Assuming a reasonable range of deposition rates for the mudstones-e.g., 1-100 m Myr -1 (Sadler et al., 1999, GeoResearch Forum.)-the studied section of the ZF was deposited millions of years after the termination of the LJE in the Onega Basin.

The section is rich in total organic carbon (TOC), with average values of 27.4±18.5 wt.% in mudstones of the OnZaP-1 and 3 cores (combined into the OnZaP section, see Methods), and 14.5±9.9 wt.% in OPH mudstones (Figure 1). These values are roughly comparable to the most organic-rich modern marine sediments (e.g., up to 21.3 wt.% on the Peru margin; Böning et al., 2004, GCA). Molybdenum, uranium (U), and rhenium (Re) concentrations in Zaonega mudstones are significantly elevated relative to other Proterozoic black shales [START_REF] Partin | Large-scale fluctuations in Precambrian atmospheric and oceanic oxygen levels from the record of U in shales[END_REF][START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF][START_REF] Reinhard | Proterozoic ocean redox and biogeochemical stasis[END_REF] , averaging 130±142 μg g -1 , 19±15 μg g -1 , and 116±84 ng g -1 , respectively, in OnZaP mudstones. In the OPH section, Mo and U average 71±92 μg g -1 and 37±50 μg g -1 , whereas overall maximum Mo, U, and Re concentrations across both sites are 1009 μg g -1 , 238 μg g -1 , and 516 ng g -1 , respectively. We focus specifically on Mo, U, and Re enrichments, as these metals have previously provided robust evidence for major shifts in Earth's redox state [START_REF] Partin | Large-scale fluctuations in Precambrian atmospheric and oceanic oxygen levels from the record of U in shales[END_REF][START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF][START_REF] Reinhard | Proterozoic ocean redox and biogeochemical stasis[END_REF] . In the OnZaP section, δ 98/95 Mo is on average 0.67±0.81‰, with a maximum of 1.49±0.14‰. Uranium isotope values in the OnZaP section range from -0.03 to 0.79‰, with an average of 0.47‰. The maximum value is the most 238 Uenriched shale value that has been reported to date (ref. 76, Yang et al., 2017, Precambrian Res.).

The primary source of these RSE to the oceans is the oxidative weathering of terrestrial RSEbearing minerals (e.g., pyrite, uraninite), and subsequent riverine transport as aqueous oxyanions [START_REF] Dunk | A reevaluation of the oceanic uranium budget for the Holocene[END_REF][START_REF] Miller | Re-assessing the surface cycling of molybdenum and rhenium[END_REF] . Hydrothermal inputs, while not fully constrained at present, are not expected to be a quantitatively significant source for any of these RSE [START_REF] Reinhard | Proterozoic ocean redox and biogeochemical stasis[END_REF][START_REF] Miller | Re-assessing the surface cycling of molybdenum and rhenium[END_REF] . The most important sink for these RSE is sequestration into marine sediments, which is influenced by the redox state of the depositional setting [START_REF] Dunk | A reevaluation of the oceanic uranium budget for the Holocene[END_REF][START_REF] Crusius | Rhenium and molybdenum enrichments in sediments as indicators of oxic, suboxic and sulfidic conditions of deposition[END_REF] . In oxic seawater, these RSE are present as recalcitrant species that tend to accumulate in the water column [START_REF] Miller | Re-assessing the surface cycling of molybdenum and rhenium[END_REF][START_REF] Anderson | Concentration, oxidation state, and particulate flux of uranium in the Black Sea[END_REF] , although Mo adsorbs to Mn(IV)-oxides under oxic conditions and can then be released in pore waters following the reductive dissolution of Mn(IV)oxides in anoxic sediments [START_REF] Algeo | Mo-total organic carbon covariation in modern anoxic marine environments: Implications for analysis of paleoredox and paleohydrographic conditions[END_REF] . In anoxic environments, these RSE are progressively converted into particle-reactive species resulting in efficient drawdown through authigenic sulphide precipitation or adsorption to organic matter [START_REF] Crusius | Rhenium and molybdenum enrichments in sediments as indicators of oxic, suboxic and sulfidic conditions of deposition[END_REF][START_REF] Anderson | Concentration, oxidation state, and particulate flux of uranium in the Black Sea[END_REF] . With regards to Mo, it has been found that most efficient Mo sequestration takes place under euxinic waters with >11 μM of sulphide (HS -) [START_REF] Helz | Mechanism of molybdenum removal from the sea and its concentration in black shales: EXAFS evidence[END_REF] . In a broad sense, the contrasting behaviour of these elements in oxic and anoxic conditions ensures that the drawdown and marine reservoir sizes are governed by global ocean redox-large reservoirs of these RSE develop in an oxic ocean, which leads to large local sedimentary RSE enrichments under anoxic conditions [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF][START_REF] Reinhard | Proterozoic ocean redox and biogeochemical stasis[END_REF] .

The U and Mo isotope systems, similar to RSE enrichments, respond to the global marine redox landscape. There is consensus that well-oxygenated oceans are characterized by heavy Mo and U isotope values (ref. 71 Using this framework, our data points to well-oxygenated oceans with high rates of primary productivity during deposition of the Zaonega Formation.

Marine RSE inventory

Enrichments in RSE (in the rest of this manuscript, RSE refers specifically to Mo, U, and Re) in the ZF far exceed anything known from pre-Neoproterozoic black shales (Figure 2). Similar values are found only in the Phanerozoic, where they have been taken as evidence for an expanded seawater RSE inventory, the direct result of pervasive ocean-atmosphere oxygenation [START_REF] Partin | Large-scale fluctuations in Precambrian atmospheric and oceanic oxygen levels from the record of U in shales[END_REF][START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF][START_REF] Reinhard | Proterozoic ocean redox and biogeochemical stasis[END_REF] . We believe that the ZF, likewise, contains evidence of a large marine RSE inventory that strongly suggests a relatively oxidized ocean-atmosphere system with a robust terrestrial oxic RSE weathering flux, in which the drawdown of RSE was limited by the relative scarcity of anoxic conditions on the seafloor. Further, U and Re enrichments-which, unlike Mo, are also drawn down in ferruginous waters [START_REF] Tribovillard | Trace metals as paleoredox and paleoproductivity proxies: An update[END_REF] -suggest that this was an episode of thorough water column oxygenation on continental shelves, rather than just a restriction of euxinic deposition, >1.4 Gyr before terminal oceanic oxygenation during the Neoproterozoic [START_REF] Och | The Neoproterozoic oxygenation event: Environmental perturbations and biogeochemical cycling[END_REF] .

While low Mo in black shales could be explained by either a small marine Mo pool, basinal restriction leading to localized Mo depletion [START_REF] Algeo | Mo-total organic carbon covariation in modern anoxic marine environments: Implications for analysis of paleoredox and paleohydrographic conditions[END_REF] , or inefficient Mo scavenging under low HS - conditions [START_REF] Helz | Mechanism of molybdenum removal from the sea and its concentration in black shales: EXAFS evidence[END_REF] , there are limited ways to explain elevated RSE concentrations in black shales of the ZF (i.e., up to several hundreds of μg g -1 of Mo throughout the >160 km-wide basin). One possibility is anomalously low sedimentation rates, which have previously been invoked to explain high TOC in U-rich black shales of the Miocene Monterey Formation, USA [START_REF] Föllmi | Phosphogenesis and organic-carbon preservation in the Miocene Monterey Formation at Naples Beach, California -The Monterey hypothesis revisited[END_REF] . Although low sedimentation rates may have played a role in concentrating TOC and RSE in the ZF, the magnitude of the RSE enrichment-consistent with Phanerozoic levels (Figure 2)-makes it unlikely sedimentation rate alone could account for it. A more plausible explanation is that the ZF was also a highly efficient RSE sink that had reliable access to a large oceanic RSE inventory [START_REF] Algeo | Mo-total organic carbon covariation in modern anoxic marine environments: Implications for analysis of paleoredox and paleohydrographic conditions[END_REF][START_REF] Sahoo | Ocean oxygenation in the wake of the Marinoan glaciation[END_REF] (see SI for additional discussion on factors controlling RSE accumulation).

The U isotope values provide an independent confirmation of highly oxygenated oceans-as noted above they are the most positive shale δ 238 U values reported to date (Figure 2), requiring both heavy seawater δ 238 U values (likely comparable to the modern value of ~-0.4‰) and a relatively full expression of the equilibrium ~1.2‰ fractionation during U reduction. Heavy seawater values only develop when oxic U burial is a major U burial flux in well-oxygenated oceans (Andersen, 2014). In addition, a significant positive U isotope fractionation (>0.6‰) strongly suggests U reduction in the water column (Andersen, 2017, Rev. Mineral. Geochem.).

Uranium reduction rates appear to scale with sulphate (or iron) reduction rates (Barnes and Cochran, 1993, GCA)-providing a link between amounts of organic matter loading in a marine system and the isotope fractionation that occurs during U burial. For reference, there is limited water column sulphate reduction in the modern Black Sea and Cariaco Basin because of low rates of productivity and maximum sediment δ 238 U values are ~0.4‰ and 0.2‰, respectively (Andersen, 2014, EPSL). Therefore, these strongly positive U isotope values, along with the RSE enrichments, are most readily explained by invoking both globally well-oxygenated oceans and markedly elevated rates of local primary productivity, the latter of which would tend to enhance the expression of isotope effects associated with U reduction.

Conflicting views on oxygenation during ZF deposition

Our interpretation that the ZF records well-oxygenated ocean conditions is ostensibly in conflict with previous interpretations derived from C, S, and Mo isotope records proposed to reflect global deoxygenation. However, much of this data can plausibly be explained through local rather than global processes. First, Kump et al. [START_REF] Kump | Isotopic evidence for massive oxidation of organic matter following the Great Oxidation Event[END_REF] postulated a post-LJE global negative δ [START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] C excursion based on depleted δ 13 C in ZF carbonates. This was interpreted to result from the widespread oxidation of organic-rich shales, which would have favoured suppressed, rather than elevated, O2 levels. However, negative δ 13 C in the ZF has subsequently been reinterpreted to reflect basinwide methanotrophy and secondary overprinting of carbonate rocks [START_REF] Qu | Methanotrophy in a Paleoproterozoic oil field ecosystem, Zaonega Formation, Karelia, Russia[END_REF] . Second, Scott et al. [START_REF] Scott | Pyrite multiple-sulfur isotope evidence for rapid expansion and contraction of the early Paleoproterozoic seawater sulfate reservoir[END_REF] argued that highly [START_REF] Dunk | A reevaluation of the oceanic uranium budget for the Holocene[END_REF] S-enriched pyrites in the upper ZF reflect a global collapse of the marine sulphate reservoir due to deoxygenation. In contrast, Paiste et al. [START_REF] Paiste | Multiple sulphur isotope records tracking basinal and global processes in the 1.98 Ga Zaonega Formation, NW Russia[END_REF] explained the same trends through basin-specific processes involving quantitative uptake of sulphate. Third, Asael et al. [START_REF] Asael | Molybdenum record from black shales indicates oscillating atmospheric oxygen levels in the early Paleoproterozoic[END_REF] inferred an ~0.7‰ seawater δ 98/95 Mo value from two other cores in the upper ZF, identical to the average in the OnZaP section. Mn-oxides, which form in oxic waters, preferentially scavenge 95 Mo and drive seawater δ 98/95 Mo to higher values. If one assumes a modern riverine input value of ~0.7‰ for the Paleoproterozoic (ref. 71), a seawater value of ~0.7‰ suggests limited Mnoxide precipitation and, thus, little O2 in the oceans. However, as outlined above, the sequestration of Mo in sediments commonly imparts a negative δ 98/95 Mo fractionation (ref. 71).

This highlights that the highest Mo isotope values of ~1.4‰ represent a minimum seawater δ 98/95 Mo estimate (see ref. 15 and SI) and that the Mo isotope data does not necessarily preclude widespread oxic conditions.

Implications of oxygenated oceans at ~2.0 Ga

RSE-replete conditions during the deposition of the ZF have important implications for Earth's global C and O cycles during the Paleoproterozoic. Tracking the evolution of Earth's redox history is also critical for understanding early eukaryote evolution. The late rise to prevalence of eukaryotes in the Neoproterozoic is thought to be linked to widespread anoxic and nutrient-poor conditions that favoured prokaryotic metabolisms [START_REF] Brocks | The rise of algae in Cryogenian oceans and the emergence of animals[END_REF][START_REF] Isson | Tracking the rise of eukaryotes to ecological dominance with zinc isotopes[END_REF] . Given the extent of ocean oxygenation and nutrient abundance inferred from the ZF [START_REF] Bekker | Oxygen overshoot and recovery during the early Paleoproterozoic[END_REF][START_REF] Lepland | Potential influence of sulphur bacteria on Palaeoproterozoic phosphogenesis[END_REF] , conditions favourable to eukaryotic diversification could have been present for much of the middle Paleoproterozoic. It is curious, then, that molecular clock analyses and microfossil evidence of eukaryote origin tend to converge after ~1900 Ma [START_REF] Parfrey | Estimating the timing of early eukaryotic diversification with multigene molecular clocks[END_REF][START_REF] Betts | Integrated genomic and fossil evidence illuminates life's early evolution and eukaryote origin[END_REF] . In the simplest sense this discrepancy supports the central importance of a unique endosymbiosis event for the emergence of eukaryotes, rather than removal or lessening of an environmental barrier.

The ZF is part of a well-recognized overall trend towards elevated RSE concentrations at ~2400-2000 Ma that has been linked to well-oxygenated oceanic-atmospheric conditions [START_REF] Partin | Large-scale fluctuations in Precambrian atmospheric and oceanic oxygen levels from the record of U in shales[END_REF][START_REF] Kipp | Selenium isotopes record extensive marine suboxia during the Great Oxidation Event[END_REF][START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF] .

However, surprisingly our ZF data extend this trend past the termination of the LJE. Furthermore, the extreme RSE enrichments hint that, instead of being in decline, atmospheric O2 abundance could still have been high perhaps several tens of Myrs after the canonical end of the LJE [START_REF] Martin | Multiple Palaeoproterozoic carbon burial episodes and excursions[END_REF] . In the modern, well-oxygenated world, the geologic O2 response time is on the order of ~2 Myrs 60 , whereas the residence times of Mo, U, and Re in the oceans are ~440, ~400, and ~780 kyrs, respectively [START_REF] Miller | Re-assessing the surface cycling of molybdenum and rhenium[END_REF][START_REF] Morford | The geochemistry of redox sensitive trace metals in sediments[END_REF][START_REF] Henderson | The U-series Toolbox for Paleoceanography[END_REF] . Therefore, if the LJE decline is linked to a decrease in O2 production as the result of diminishing Corg burial by the end of the LJE, then O2 would be expected to be significantly attenuated by this time. Instead, our observations provide evidence for the opposite.

Our observations provide empirical evidence for models (e.g., ref. 18) that decouple the strongly positive carbonate carbon isotope values of the LJE from enhanced organic carbon burial. Further, if highly oxidized conditions continued up to ~2.0 Ga-or, alternatively, if the Paleoproterozoic experienced large fluctuations in O2 conditions, the ZF being a zenith-it follows that O2 abundance was governed by more complex controls than a simple Corg burial mass balance, as outlined in Bekker and Holland [START_REF] Bekker | Oxygen overshoot and recovery during the early Paleoproterozoic[END_REF] . This highlights the importance of exploring alternative mechanisms that could account for Paleoproterozoic O2 accumulation and the need for refined chemostratigraphic and geochronological studies focused on this key interval of Earth's history.

Methods

Materials

The material for this study comes from drill cores in the Onega Basin that intersect the ZF-a relatively well-preserved 1500 m-thick succession of organic-rich mudstones and carbonates, interlayered and intersected with igneous units including lavas, tuffs, and sills (see SI for a detailed geological setting) [START_REF] Melezhik | Karelian shungite-an indication of 2.0-Ga-old metamorphosed oil-shale and generation of petroleum: geology, lithology and geochemistry[END_REF][START_REF] Martin | Multiple Palaeoproterozoic carbon burial episodes and excursions[END_REF][START_REF] Melezhik | A giant Palaeoproterozoic deposit of shungite in NW Russia: genesis and practical applications[END_REF][START_REF] Melezhik | The Onega Basin[END_REF] . The 60 m long cores OnZaP-1 and OnZaP-3 were drilled 500 m apart in Shunga village, close to drill core FAR-DEEP 13A [START_REF] Paiste | Multiple sulphur isotope records tracking basinal and global processes in the 1.98 Ga Zaonega Formation, NW Russia[END_REF][START_REF] Črne | Zaonega formation: FAR-DEEP hole 13A. in Reading the Archive of Earth's Oxygenation[END_REF] . The 3500 m long OPH core was drilled ~70 km to the south. Paiste et al. [START_REF] Paiste | Multiple sulphur isotope records tracking basinal and global processes in the 1.98 Ga Zaonega Formation, NW Russia[END_REF][START_REF] Paiste | Reconstructing the Paleoproterozoic sulfur cycle: Insights from the multiple sulfur isotope record of the Zaonega Formation[END_REF] correlated the three cores based on C isotopes, trace metal enrichments, P concentrations, and a distinct massive P-rich dolomite unit that occurs throughout the Onega basin [START_REF] Črne | Petrography and geochemistry of carbonate rocks of the Paleoproterozoic Zaonega Formation, Russia: Documentation of 13C-depleted non-primary calcite[END_REF] .

The lithology of this section is characterized by alternating dolomite to calcite-rich carbonates and exceptionally organic-rich mudstones that are intersected by silica or pyrobitumen veins (Figure 1). In the OnZaP section, the interval from the bottom of the core to 53 m depth is dominated by highly organic-rich mudstones with relatively few carbonate beds; the 53-33 m interval is mostly dolostone and contains a distinctive dolomite unit; and the 33-1.7 m interval consists of grey mudstones and marly carbonate beds [START_REF] Paiste | Multiple sulphur isotope records tracking basinal and global processes in the 1.98 Ga Zaonega Formation, NW Russia[END_REF] . Two sets of samples from the OnZaP section were analyzed in this study. Set MSP0001 consists of 135 samples that were taken at roughly 1 m intervals and is identical to that used in Paiste et al. [START_REF] Paiste | Multiple sulphur isotope records tracking basinal and global processes in the 1.98 Ga Zaonega Formation, NW Russia[END_REF] . Set MSP0010 contains 79 samples that more specifically targeted RSE and organic-rich intervals. A set of 89 samples were analyzed from the 1060-1230 m interval of the OPH core that is roughly equivalent to the OnZaP section.

Elemental concentrations

Total organic carbon content of the MSP0001 and OPH samples were adapted from Paiste et al. [START_REF] Paiste | Multiple sulphur isotope records tracking basinal and global processes in the 1.98 Ga Zaonega Formation, NW Russia[END_REF] and were measured from powdered aliquots using a LECO SC-444 analyzer at the Geological Survey of Norway (NGU), Trondheim, Norway. Detection limit was 0.1 wt.% and precision better than 10%. For the MSP0010 set, TOC was measured at the Pôle Spectrométrie Océan, European Institute for Marine Studies (IUEM), Brest, France. Dried and powdered samples were combusted in ceramic beakers at 500 °C for 24h and loss of mass on ignition was determined. Repeat measurements of 7 samples generally differed by <1 wt.%.

Major and trace element compositions of MSP0001 and OPH samples were determined at Acme Labs, Bureau Veritas Commodities Canada Ltd. The samples were first pulverized and combusted to remove organic carbon. For major elements, the sample was fused in LiBO2 flux, then digested, whereas minor elements were analyzed from full digests in HNO3, HClO4, and HF.

Measurements utilized inductively coupled plasma optical emission spectrometry (ICP-OES) or mass spectrometry (ICP-MS). Average relative standard deviation was less than 5% for all elements. MSP0010 element composition was determined at IUEM. Samples were pulverized in a tungsten carbide crusher and an agate mill, then combusted at 500 °C for 24h to remove organic carbon. Major element concentrations were measured on a Jobin Yvon Horiba Ultima 2 ICP-OES after digestion overnight at 80 °C in concentrated HF and HNO3 and neutralization with boric acid to retain Si. For trace element concentrations, digestion was performed in a class 1000 clean laboratory using distilled acids-samples mixed with concentrated HF and HNO3 were heated to 80 °C overnight, then allowed to evaporate; this digestion was then repeated with concentrated aqua regia, after which the sample was taken up in 6M HCl. Trace element concentrations were [START_REF] Kendall | The stable isotope geochemistry of molybdenum[END_REF][START_REF] Mclennan | Relationships between the trace element composition of sedimentary rocks and upper continental crust: Trace element composition and upper continental crust[END_REF][START_REF] Tissot | Uranium isotopic compositions of the crust and ocean: Age corrections, U budget and global extent of modern anoxia[END_REF] and the dotted horizontal line is the phosphorous-rich mudstone-dolostone contact used for intra-basinal correlation (see Methods and SI). Boxplots below the Mo and U isotope plots represent the distribution of previously published ZF data from lower in the succession [START_REF] Asael | Molybdenum record from black shales indicates oscillating atmospheric oxygen levels in the early Paleoproterozoic[END_REF][START_REF] Asael | Coupled molybdenum, iron and uranium stable isotopes as oceanic paleoredox proxies during the Paleoproterozoic Shunga Event[END_REF] . Gray horizontal bands are partly silicified and calcified intervals due to fluid alteration [START_REF] Paiste | Multiple sulphur isotope records tracking basinal and global processes in the 1.98 Ga Zaonega Formation, NW Russia[END_REF] . While these imply the presence of secondary fluids that could have mobilized RSE, microscale RSE distribution confirms the primary nature of the RSE enrichments (see SI).

Figure 1 :

 1 Figure 1: Lithology and geochemistry of the Zaonega Formation. Total organic carbon (TOC) and redox sensitive metal (RSE) profiles, and Mo and U isotope ratios (δ 98/95 Mo and δ 238 U) for the OnZaP and OPH sections. Vertical dashed lines are average crustal values 71-73 and the dotted

Figure 2 : 20 Figure 1

 2201 Figure 2: Secular trends in redox sensitive element concentrations from anoxic shales: Zaonega Formation data plotted on compilations from literature. (a) Changes in δ 13 Ccarb ratio through time, modified from Karhu 74 . (b) Molybdenum concentrations from Robbins et al. 75 . (c) Uranium concentrations from Partin et al. 11 . (d) Uranium isotope ratios from Wang et al. 76 , Yang et al. 2017. (e) Rhenium concentrations from Sheen et al. 13 .

  and Andersen et al., 2017, Rev. Mineral. Geochem.). It is, however, difficult to gauge when a sedimentary archive captures seawater U and Mo isotope values. For instance, it is easiest to capture seawater Mo isotope values in shales deposited in isolated basins with high HS -levels (ref. 50), but euxinic marine basins often deviate from the seawater isotope composition (ref. 71)-δ 98/95 Mo in euxinic shales is therefore often interpreted as a minimum estimate for seawater (ref. 15). Uranium isotope values in shales will be closest to seawater values in oxic sediments and up to ~1‰ heavier than seawater in anoxic and high productivity settings (Andersen et al., 2014, EPSL) (see SI for detailed discussion on Mo and U isotope values).

Acknowledgements E. Ponzevera, A. De Prunelé, M.L. Rouget, and K. von Gunten are thanked for help with trace element analyses. KM, KK, AL, TK, PP and KP were supported by the Estonian Science Agency grant PRG447. KM was further supported by the Ministry of Education and Research of Estonia mobility grant within Archimedes Foundation's The Kristjan Jaak Scholarship program "Doctoral Study Abroad"; KP was supported by the Research Council of Norway through its Centres of Excellence funding scheme grant No. 223259; KOK was supported by the Natural Sciences and Engineering Research Council of Canada Discovery grant (RGPIN-165831).

Data availability

The novel ZF geochemical data presented here are included in the supplementary files.

measured on a Thermo Scientific Element 2 ICP-MS calibrated against commercial multi-element standards and digested geostandards (e.g. BHVO-2) were Samples with a ratio of (CaO+MgO)/(SiO2+Al2O3) > 0.2-approximating >20 wt.% carbonate content-were excluded when calculating Mo, U, and Re averages, since compilations of RSE in black shales generally only include shale samples [START_REF] Partin | Large-scale fluctuations in Precambrian atmospheric and oceanic oxygen levels from the record of U in shales[END_REF][START_REF] Sheen | A model for the oceanic mass balance of rhenium and implications for the extent of Proterozoic ocean anoxia[END_REF][START_REF] Reinhard | Proterozoic ocean redox and biogeochemical stasis[END_REF] . An exception was made for samples with >5 wt.% TOC, since such sediments would have played a role in trace metal cycling regardless of their mineralogical composition.

Mo isotopes

Mo isotopes were measured from fully digested OnZaP MSP0010 samples that were purified via column chromatography according to ref. 68. A 97 Mo-100 Mo double spike was employed [START_REF] Siebert | Determination of molybdenum isotope fractionation by double-spike multicollector inductively coupled plasma mass spectrometry[END_REF] and the isotopes were measured on a Thermo Scientific Neptune multi-collector ICP-MS at IFREMER, Brest, France. Data is expressed relative to NIST SRM 3134 = 0.25‰ [START_REF] Nägler | Proposal for an international molybdenum isotope measurement standard and data representation[END_REF] . Detailed methodology is provided in the SI.

U Isotopes

Sample for U isotopes were sequentially digested in a mixture of 3 mL HNO3 and 1 mL HF at 100ºC for 24 hours, then aqua regia at 95ºC for 24 hours. Following each digestion step the sample was evaporated to dryness. Sample resides were taken up in 5 mL of 3 M HNO3 at 70ºC.

All sample preparation was performed in a

Competing interests

The authors declare no competing interests.


Tat Loon Chng

Maya Naphade

Benjamin Goldberg

Igor Adamovich

Svetlana M Starikovskaia

Towards Ultrashort Pulse Generation from Nanosecond Laser Sources for Electric Field Measurements

Introduction

A laser diagnostic based upon the phenomenon of second harmonic generation has been re-developed recently for the purpose of electric field measurements in gas plasmas [1][2]. More commonly known as electric field induced second harmonic generation (E-FISH) [3], the relatively straightforward implementation of this technique has begun to find application in real-life plasma actuators [4][5]. Since the signal generation is intrinsically intensity-driven, this naturally favours short laser pulses on the order of hundreds of picoseconds (ps) or less. The use of high-energy, nanosecond pulses (ns) is in principle possible, but increasing the energy indefinitely is constrained by the occurrence of laser-induced breakdown. Furthermore, ns pulses are limited in time-resolution, which may prove important when diagnosing fields associated with ns type discharges. Unfortunately, ultrashort pulsed laser systems are not easily accessible and costs can often be prohibitive.

In this work, we attempt to use a Pockels cell to 'slice' a ubiquitous 1064 nm YAG laser beam into a pulse of shorter duration, with a view to using this shorter pulse for useful E-FISH measurements in ns discharges.

Experimental setup

Fig. 1: E-FISH setup

The experimental setup consists of 3 main parts. The high voltage pulses and electrode system used for demonstrating the E-FISH measurements, the optical layout for the pulse-slicing, and that for the E-FISH measurements. An FID model number FPG 12-1NM high-voltage (HV) generator delivers positive polarity, 9.4 kV amplitude voltage waveforms with a rise time of ~8 ns and a full width at half maximum (FWHM) of ~30 ns at a repetition rate of 10 Hz. This voltage is applied across two rectangular, flat plate electrodes separated by a gap of 4 mm as shown in Fig. 1.

The Pockels cell used in these experiments was supplied by Leysop Ltd and consists of a UV BBO crystal housed in a custom-made casing with electrical leads. Applying a HV pulse to the Pockels cell alters the birefringence of the BBO crystal for the duration of this pulse, effectively rendering the cell as an 'electronic waveplate'. By placing the Pockels cell between two crossed polarizers and timing it such that it is triggered in coincidence with a laser pulse, a shorter laser pulse may be obtained (see Fig. 2) [6]. The energy of the sliced pulse is strongly dependent on the polarization rotation that can be produced by the Pockels cell. As seen in the top schematic of Fig. 2, for a fully transmitting cell, performance is best when half-wave rotation can be effected. Unfortunately, the present cell was designed for 205 nm light and limits the amount of polarization rotation that can be obtained at 1064 nm. Furthermore, the possibility of crystal damage at this off-optimum wavelength also restricts the energy input into the cell. To mitigate this problem, a double-pass arrangement is implemented, results of which are shown in the following section.

The E-FISH setup given in Fig. 1 is essentially similar in principle to that described in [2]. A Quanta-Ray Lab-230 YAG laser system is used as a source of 1064 nm, 50 mJ, ~15 ns pulses at repetition rate of 10 Hz. These pulses are reduced to about 2.8 ns and 2 mJ after passing through the double-pass pulse slicing layout given in the bottom panel of Fig. 2. This beam is focused with a 30 cm lens into the center of the electrode gap and collimated with a 25 cm lens. A dichroic mirror and dispersive prism separate the residual 1064 nm beam and the 532 nm second harmonic (SHG) signal. This 532 nm light is then focused with a 10 cm lens into a photomultiplier with an attached 532 nm bandpass filter. Part of the residual 1064 nm beam is picked off and directed onto a photodiode as a way of monitoring the laser intensity.

The intensity of the 532 nm, E-FISH signal, 𝐼𝐼 𝑦𝑦 (2𝜔𝜔) , is

given by,

𝐼𝐼 𝑦𝑦 (2𝜔𝜔) ∝ �𝜒𝜒 𝑦𝑦𝑦𝑦 (3) (2𝜔𝜔, 0, 𝜔𝜔, 𝜔𝜔)𝑁𝑁𝐸𝐸 𝑦𝑦 𝐼𝐼 𝑦𝑦 (𝜔𝜔) 𝐿𝐿� 2 �𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 � 𝐿𝐿 𝐿𝐿 𝑐𝑐 �� 2 (1)
In Eq. ( 1), 𝐸𝐸 𝑦𝑦 is the vertical (y-axis) component of the electric field, 𝐼𝐼 𝑦𝑦 (𝜔𝜔) is the intensity of the vertically polarized 1064 nm laser beam, 𝜒𝜒 𝑦𝑦𝑦𝑦 (3) is the third-order nonlinear susceptibility, N is the gas number density, L is the confocal beam parameter, 𝐿𝐿 𝑐𝑐 is the coherence length where

𝐿𝐿 𝑐𝑐 = 2 𝛥𝛥𝛥𝛥
and Δk is the difference between the fundamental and the second harmonic wave vectors.

Results and discussion

Fig. 3 shows the results obtained with both the single and double-pass arrangements. The double-pass option produces about 50 times more energy without any detriment to the pulse width. This encouraging result supported the use of this multi-pass configuration for all the results that follow. The overall energy conversion efficiency in this case is estimated to be about 4%. To evaluate the viability of these shorter pulses, the E-FISH signal was measured at various time instants in a Laplacian field produced by the HV ns waveforms. This was achieved by stepping through the delay of the main laser beam relative to the applied HV pulses. Fig. 4 shows the voltage evolution obtained using the E-FISH technique with the sliced pulses. Measured signals are calibrated against corresponding values measured by a custom-made back current shunt (BCS) installed in the cable connected to the high-voltage electrode. It is clearly seen that these shorter pulses are capable of generating signals which follow the fast rise and fall of the voltage waveforms.

Conclusion

The use of an off-optimum design Pockels cell to generate short pulses for electric field measurements has been successfully demonstrated. Significantly shorter and higher energy pulses can be readily achieved with a welldesigned, commercial package and at a fraction of the cost of a typical ultrashort pulse laser system. This economical solution may pave the way for performing useful measurements for a variety of laser-based diagnostics.

Fig. 2

 2 Fig. 2 Pulse slicing optical schematic. Top: Single-pass, Bottom: double-pass.

Fig. 3

 3 Fig. 3 1064nm laser pulse waveforms obtained with (a) single-pass and (b) double-pass configurations.

Fig. 4

 4 Fig. 4 Laser intensity normalized E-FISH signals plotted against back current shunt data. Dividing these voltages by the interelectrode gap distance gives the electric field strength.

Acknowledgements

This work was partially supported by the French National Research Agency, ANR (ASPEN Project), LabEx Plas@Par and the French-Russian international laboratory LIA KaPPA 'Kinetics and Physics of Pulsed Plasmas and their Afterglow'. The support of Prof. Adamovich by the Ecole Polytechnique Gaspard Monge Visiting Professor (GMVP) Program is gratefully acknowledged.