Kenichi Koizumi

Makoto Hatakeyama

Mauro Boero

Katsuyuki Nobusada
email: nobusada@ims.ac.jp

Hirokazu Hori

Taku Misonou

Shinichiro Nakamura

How seaweeds release the excess energy from sunlight to surrounding sea water

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

Introduction

Sunlight is the primary energy source of living organisms on the earth. and, in this context, photosynthesis is one of the most important metabolism processes. It consists in the synthesis of nutrients from the environmental carbon dioxide and water, catalyzed by visible sunlight. Yet, the solar spectrum is wider than the simple visible light and includes also high-energy ultra-violet (UV) radiation (UVR), harmful to any living organism. Indeed, UVR can damage to various extents biomolecules and, specifically, nucleic acids (DNA and RNA) in the cell through direct and indirect pathways. To cope with this problem, several evolution strategies have been adopted by nature to protect life from the UVR damage. A clear understanding of these strategies is now becoming increasingly important also for the sunlight-exposed biosphere because of the constant increasing of UVR on Earth, particularly enhanced in the last decades, resulting from the disruption of the ozone layer in the stratosphere mainly due to an anthropogenic factor.

The development and accumulation of UV absorbing substances inside the cells, such as mycosporine-like amino acids (MAA), can be regarded as a powerful physiological solution and a fundamental step against UVR adopted by nature in general and algae in particular. [1][2][3] MAA was found in a wealth of organisms, particularly n green-blue 4,5 and red algae, 6 living in environment highly irradiated by sunlight and display a maximum absorption in the wavelength range between 300 and 350 nm, thus absorbing mainly UV-A and UV-B. 7 In fact, MAA is ubiquitous in various aquatic and terrestrial vegetal organisms, namely cyanobacteria, algae and fungi; moreover it is also present in some aquatic animals. 8 The function exerted by MMA extends beyond that of an UV protector, and includes also other physiological roles as an antioxidant, an osmolyte and as an agent preventing desiccation and thermal stress. [9][10][11] This rich variety of functions make MMA a fundamental compound for life in general and its study a fundamental multidisciplinary target from photocatalysis to advanced biomaterials. Concerning its primary function, the UV absorption characteristics of MAA have been extensively studied with respect to their physiological and ecological importance. Instead, research activities focused more on the physico-chemical aspects and energy transfer mechanisms of MAA molecules are rare and still at a pioneering stage. Yet, the efficiency of MMA in energy adsorption upon UVR is remarkable and worthy of the attention of research efforts at the crossroad among physics, chemistry and biology. In fact, the molar extinction coefficient of MMA can reach extraordinary values up to a range 28000 -50000 M -1 cm -1 with unsurpassed sunscreen properties. 12,13 Among the studies available to date, Misonou and coworkers [START_REF] Misonou | Proceedings of 3rd International Scientific-practical conference on Marine coastal ecosystem[END_REF] have shown that MAA is able to quench the excited state of a pyrimidine nucleotide, a typical component of nucleic acids. While MAA is known to possess a high stability even under relatively strong UV irradiation, [START_REF] Whitehead | [END_REF] the absorbed energy is expected to be dissipated very rapidly in the surrounding water environment, [START_REF] Misonou | Proceedings of 3rd International Scientific-practical conference on Marine coastal ecosystem[END_REF][16][17][18] although the mechanism responsible for this dissipation is largely unknown. Furthermore, MAA has a low level of fluorescence (1 × 10 -4 ), [16][17][18][19] a fact that makes hard unraveling the mechanism via experimental probes.

The scope of the present study is then to address the fundamental question of the excess energy dissipation and to shed some light into the basic mechanism responsible for this process. To this ometries are the most stable. Energetically they are rather close to each other and correspond to the configurations labeled as 2 and 8 in Figure 1. Details will be discussed in the next paragraph. After the determination of these two ground state stable configurations, we computed the corresponding excited-state geometries via a single-excitation configuration interaction (CIS) approach. In all these quantum chemical calculations, a localized 6-31G(p,d) basis sets 41 was adopted and the effect of the solvating water was kept into account by the inclusion of a polarizable continuum model (PCM). 42,43 For this set of static calculations, the GAUSSIAN09 package [START_REF] Frisch | Revision E.01[END_REF] was used. All subsequent dynamical simulations were conducted at the FPMD level within the Car-Parrinello molecular dynamics (CPMD) [START_REF] Car | [END_REF] method. Valence electrons wavefunctions were expanded on a plane-wave basis set with an energy cut-off of 80 Ry and the sampling of the Brillouin zone was restricted to the Γ point. All core-valence interactions were described by analytical norm-conserving pseudopotentials (PPs) of the Hartwigsen-Goedecker-Hutter (HGH) type. [47][48][49] The generalized gradient approximation (GGA) adopted in all FPMD simulations is the one composed of the Becke exchange 39 and the Lee-Yang-Parr-correlation 40 (BLYP). The Grimme's empirical formula 50 was included to account for van der Waals interactions, generally missed by any GGA. The efficiency of such a computational set-up has been extensively benchmarked over the years on a wealth of aqueous systems [51][52][53][54] and can be regarded as sufficiently accurate.

The porphyra-334 molecules that were found as stable, namely configurations 2 and 8, were solvated in a periodically repeated supercell of liquid water at standard density (1.0 gcm -3 ) and temperature (300 K). Given the encumbrance of these porphyra-334 structures, 120 and 125 water molecules were necessary to fill a cubic simulation cell of lateral size of 17.5 Å . The construction of these systems started by a classical molecular dynamics equilibration of the solvent water alone. Then, to insert the porphyra-334 molecular structures, a sufficient space was cleaved by eliminating all H 2 O molecules falling inside the van der Waals radii of the optimized porphyra-334 geometries in their excited states (see figure S1 in the SI). Then the system was re-equilibrated for 10 ps within a CPMD approach at 300 K in the canonical NVT ensemble with a Nosé-Hoover temperature control scheme. [55][56][57] Initially, the structure of the porphyra-334 molecule was kept fixed to its excited state configuration, whereas the surrounding water was equilibrating. After the equilibration of the solvent, all constraints on the porphyra-334 geometry were released and the system was allowed to evolve dynamically on the underlying potential energy surface (PES). At this point, in order to eliminate any thermodynamical bias, the thermostats were suppressed and all simulations continued in the microcanonical NVE ensemble. Such a set-up accounts for the initial excitation of the solute upon sunlight adsorption, a process which is generally rapid enough to have little or no immediate influence on the surrounding solvent. Then, to reproduce the deactivation processes, initial velocities are attributed to the atoms of the porphyra-334 molecules to follow the non-radiative transition process that NVE simulations mimic, schematically illustrated in Figure 2.

The excited porphyra-334 molecule is expected to relax to the ground electronic state in a time much shorter than the accompanying atomic relaxation, 58,59 thus all the electron excitation energy is converted into an initial excess kinetic energies accumulated in the porphyra-334 molecule.

Since CPMD simulations are done on the ground state electronic configuration, to reproduce such an initial excess kinetic energy, random velocities, sampled from a Maxwell-Boltzmann distribution corresponding to a typical 334 nm photon wavelength excitation, are attributed to the atoms of the porphyra-334 molecule to mimic these conditions. This amounts to an initial excess (kinetic) energy of K = 3.7 eV redistributed among the porphyra-334 atomic constituents and, hence, resulting in a strong heating of about 624 K according to the standard formula K = 3/2Nk B T . Here, k B is clearly the Boltzmann constant and T , N and K indicate the temperature, the total number of atoms, and the total kinetic energy, respectively. Our CPMD simulations, then, account for the relaxation processes involved in the release of this extra kinetic energies of the porphyra-33 to the surrounding solvent water.

After an initial 10 ps equilibration of the solvent, with the porphyra-334 kept fixed, we continued the CPMD dynamical simulation observing a continuous breaking and formation of H-bonds around the solute molecule. This implies that the H-bond network around the porphyra-334 is highly dynamical and can rearrange on a ps time scale. In order to account for these different H-bond configurations and to increase the statistics, for both configurations 2 and 8 we selected at intervals of 10 ps seven different configurations for each one. These time intervals are long enough to ensure that the sampled configurations are uncorrelated and the net result was a set of fourteen different configurations (seven for configuration 2 and seven for configuration 8) from which independent simulations of the porphyra-334 relaxation were started. To this aim, for each simulation a randomization of the initial velocities was done according to the procedure described above. Each one of these fourteen independent simulations lasted for 9 ps, a time that, although short, turned out to be sufficient to observe the total release of the excess kinetic energy initially stored in the porphyra-334. The process was then monitored and analyzed in term of mass-weighted autocorrelation functions for which details are given in the following subsection.

Analysis of the Simulations Data

To analyze the relaxation processes in which an energy transfer occurs from the porphyra-334 molecule to the surrounding solvating water, we resorted to the methodology proposed by Yagasaki and Saito. [60][61][62] Such a procedure allows for a detailed inspection of the evolution in time of the kinetic energy transfers from the heated porphyra-334 to the solvent. Specifically, the extra kinetic energies accounting for the deexcitation of the porphyra-334 is partitioned into separate contributions, namely the intermolecular motion, the internal rotation, and the vibrational modes.

Explicitly, the excess kinetic energy is written as

∆C(t , T ) = 1 2 ∑ i m i v v v i (t + T )v v v i (T ) NE - 1 2 ∑ i m i v v v i (t + T )v v v i (T ) EQ , (1)
The superscript "NE" (non-equilibrium) indicates the ongoing relaxation process in which extra kinetic energy of the porphyra-334 molecule is gradually transferred to surrounding H 2 O molecules. The superscript "EQ", instead, refers to the dynamics in which both the porphyra-334 and the solvating water are at the equilibrium, i. e. the system is an aqueous solution at 300 K in the present case. The mass and the velocity vector of the ith atom are indicated as m i and v v v i , respectively. T is the selected waiting time starting from which the calculation is done and t is the running variable indicating the time evolution from T onward. Hence, the first and second terms in the left-hand part of equation ( 1) are not directly kinetic energies (although resembling) but mass-weighted auto-correlation functions. A real Fourier transform, i. e. making use of the cosine only, gives

∆C(ω t , T ) = 1 π ∞ 0 dt cos(ω t t ) ∑ i m i v v v i (t + T )v v v i (T ) NE - 1 π ∞ 0 dt cos(ω t t ) ∑ i m i v v v i (t + T )v v v i (T ) EQ .
(

) 2
This equation can be interpreted as the frequency dependent difference between the kinetic energies of the non-equilibrated dynamics of the solvated excited porphyra-334 and the one of the equilibrated system. The above difference of frequency dependent kinetic energies can then be splitted into two separate contributions,

∆C(ω t , T ) = ∆C tra (ω t , T ) + ∆C rot (ω t , T ), (3)
where ∆C tra is the part accounting for the translational motion of the solvating H 2 O molecules, thus can be directly computed from the kinetic energies of the center of mass of the water molecules only. ∆C rot indicates, instead, the excess kinetic energy due to rotational and vibrational motions and can be calculated by subtracting the contributions. of the kinetic energies of the water molecules from the total kinetic energy of the whole solvent. On these grounds, we can define the partial excess kinetic energies (PEKEs) distributed by the five inter or intra molecular motion modes as

∆K HR (T ) = ω vib ω rot dω t ∆C rot (ω t , T ), (4a)
∆K LR (T ) = ω rot 0 dω t ∆C rot (ω t , T ), (4b)
∆K V IB (T ) = ∞ ω vib dω t ∆C rot (ω t , T ), (4c)
∆K HT (T ) = ∞ ω tra dω t ∆C tra (ω t , T ), (4d)
∆K LT (T ) = ω tra 0 dω t ∆C tra (ω t , T ). (4e)
Here, ω rot and ω tra are set to 600 and 110 cm -1 , respectively, according to the prescription proposed by Yagasaki and Saito. [60][61][62] Analogously, ω vib is set to 1250 cm -1 , a value able to fully separate the intramolecular vibrational contributions from the rotational ones. The type of motion represented by HR, LR, HT, and LT are sketched for clarity in Figure 3. With this algorithm, the excess kinetic energy is separated into five contributions (the fifth being the vibrational motion contribution) and these can discriminate between inter and intra molecular motions. From a strictly physical point of view, we have to remark that these PEKEs do not have the dimensions of a kinetic energy, but rather those of a kinetic energy (kcal/mol) multiplied by a wave number (cm -1 ).

Results and Discussions

Structural optimizations of the porphyra-334 in its ground and excited states To prepare the system for subsequent dynamical simulations, geometry optimizations of the porphyra-334 molecule in its ground state were done within the DFT framework, by using a B3LYP functional to treat the exchange and correlation interactions. By assuming that the porphyra-334 is in a neutral charge state, the effects of a surrounding solvent were accounted for by a PCM model 42,43 in all structural relaxations. As mentioned in the former paragraph, to account for all possible protonation states of hydroxy groups and nitrogen atoms of the porphyra-334, nine different structures were considered (see Figure 1). Our search for stable configurations has shown that the configurations labeled as 2 and 8 in Figure 1 are characterized by an energy difference of 0.15 eV, the lower reference value being the configuration 8, and they are the most stable among all the possibilities analyzed. In all the cases in which an NH 2 group is present in the molecule and when carboxyl ions are located in proximity of these NH 2 groups, one of the two H of the NH 2 moiety departs and migrates toward the carboxyl group nearby. The net result is that the geometries of configurations 3 and 4 in Figure 1 transform spontaneously into configurations 2 and 8, confirming again their better overall stability. Energetically speaking, above the two configurations 2 and 8, we could locate the structure indicated as configuration 1. This specific structure is higher in energy by about 0.45 eV with respect to configuration 8 and 0.30 eV with respect to configuration 2, the most stable being the former. Hence, on these grounds, we selected configurations 2 and 8 as possible porphyra-334 structures for subsequent excited state calculations and CPMD simulations.

The optimized geometries of these two configurations are sketched in Figure 4.

Starting from the ground state optimized geometries for configurations 2 and 8, we re-optimized these two structures in their respective excited states at the CIS level by including a polarizable continuum model 42,43 to account for solvent effects. In both cases singlet states were considered by imposing a total multiplicity 2s + 1 = 1. The results in terms of structural changes between the excited and the ground state geometries are shown in Figure 5 ( The geometries of the excited states can be regarded as the ones corresponding to the point indicated as B in Figure 2 and they are starting structures for all subsequent simulations of the relaxation processes. A comparison of the excited and ground state geometries shows that the cyclohexane moieties become slightly bended upon electron excitation. Analogously, all carbon backbones around hydroxy and methyl groups show noticeable distortions with respect to the ground state configurations. The root mean square deviations (RSMDs) between the ground and the excited states turn out to be 0.51 Å and 0.52 Å for the configurations 2 and 8, respectively.

First principles molecular dynamics simulations

The excited state geometries of configurations 2 and 8 as provided by the CIS calculations were used as starting conformations for CPMD simulations. As explained in the section dedicated to the computational details, after solvating these structures in liquid water, we proceeded to a reequilibration of the solvent keeping the porphyra-334 fixed to the excited geometry. Once this equilibrium was attained, we attributed random velocities to the atoms of the porphyra-334 molecule sampled from a Maxwell-Boltzmann distribution corresponding to a temperature of 624 K. This specific temperature, as explained above, accounts for the energy gained by the solute upon photon adsorption. With the two systems prepared in such a way, all constraints were released and dynamical simulations performed in the microcanonical NVE ensembles. According to the analysis described by equation ( 2), by setting the waiting time T to zero when the equilibrium is reached we calculated the frequency dependent kinetic energy of the solvating water molecules, specifically the real Fourier transform of the mass-weighted velocity autocorrelation function. The result of this analysis is summarized in Figure 6. The disentangling of the frequency dependent kinetic energy from the contributions due to molecular translations, rotations and vibrations was done according to the definitions given in the computational methods section. This analysis shows that the frequency dependent kinetic energy of the translational motion has a pronounced peak centered at 100 cm -1 accompanied by a shoulder around 200 cm -1 (green line in Figure 6). The rotational contributions, on the other hand, are characterized by a peak centered at 600 cm -1 . Peaks at 1600 cm -1 and 3200 cm -1 represent, as usual, the bending and stretching molecular vibrations of liquid water and are clearly well separated by the selected threshold ω vib = 1250 cm -1 discussed in the previous paragraph. These spectral distributions reproduce with appreciable accuracy the vibrational density of states provided by both classical simulations and experiments 60 and represent a benchmark for our DFT-based CPMD simulations.

Frequency (cm -1 ) kcal/mol stretching the preferred channel for the dissipation of the excess kinetic energy of the porphyra-334 molecule. In fact, between the two configurations 2 and 8, the increase of the HT PEKE contributions as a function of time is remarkably larger in the latter case. Only minor or nearly negligible contributions can be appreciated by observing the trend of the PEKEs due to the LT, LR, HR, and vibrational contributions. It can be inferred that the rather large kinetic energy gained upon photoexcitation directly activates the higher translational modes, hence HT rather than the lower LT mode Most likely this is due to the fact that the HT mode can accept a larger amount of excess energy with respect to the other modes, namely LT, LR, and HR.

Since the HT motion involves the collective shrinking and elongation of intermolecular Hbonds, this dynamical change of the distance between H 2 O molecules can result in the breaking (and reforming) the local H-bond network which, in turn, stabilize the hydration shell of the porphyra-334 molecule. More precisely, a decrease in the number of H-bonds destabilizes temporarily the system and produces an increase in the potential energy of the solvent water. To make this statement more quantitative, we analyzed the evolution of the total number of H-bonds during all the 9 ps dynamical simulations, considering both the H-bond network of the solvent and the Hbonds of the solute with the surrounding water. To this aim, the H-bonds between H 2 O molecules and between the solvent water and the porphyra-334 solute were monitored by using the Visual Molecular Dynamics (VMD) tool. 63 A snapshot of the continuous H-bond network of the system is shown in the SI (Figure S1) and its local evolution along the dynamics is reported in panels (a)

and (b) of Figure 9. A simple visual inspection of Figure 9 shows that the number of H-bond between solvent and solute is highly variable as the dynamics proceeds and this, in turn, influence the whole solvent.

Such a process has been quantified by monitoring H-bond formation and disruption. For the solvent only, the result of this analysis is summarized in Figures (10)) for the systems whose simulations were started from configurations 2 and 8, respectively. This result represents the average over all the seven independent simulations. Details for each single simulations are reported in the SI, Figures S4 andS5.

By looking at the trend of the rolling average, in the case of simulations started from the prophyra-334 in the configuration 2 we can observe a gradual decreases in the number of H-bonds from an initial value of ∼75 to a final value of ∼60. An analogous trend characterized also the evolution for the system started from configuration 8, where, again, a drop in the number of Hbonds can be observed as the simulation progresses. The overall conclusion that can be drawn is that the activation of the intermolecular HT motion destabilized and disrupts to some extent the

Conclusions

After determining the possible stable configurations of a porphyra-334 molecule in its ground state, we inspected the structural changes occurring upon photon adsorption via precise quantum chemi-20 cal calculations at the CIS level. On these grounds, first principles molecular dynamics simulations were used to inspect in detail the release of the excess kinetic energy accumulated in the excited porphyra-334 toward the surrounding environment represented by water. Our simulations were able to reproduce the nonradiative transition processes involved in such an energy release, clarifying both the mechanism and the preferred channel along which such an excess energy is released when the porphyra-334 relaxes back to its ground state. The disentangling of the different contribution in terms of partial excess kinetic energies attributed to the various collective modes of the solvent has shown that the main pathway along which the kinetic energy transfer occurs is the hydrogen bonds stretching (HT mode) among water molecules. Indeed, this is the mode that induces a minimal perturbation and distortion of the continuous H-bond network, as opposed to molecular rotations (HR, LR) and bending (LT). These HT stretching mode can result in instantaneous breaking, but an equally instantaneous reformation of existing H-bonds within the continuous network of the tetrahedrally coordinated H 2 O constituents of the solvent. A precise monitoring of the number of H-bonds during the dynamical evolution of the system toward the equilibrium has shown that this number decreases gradually by ∼30 in the solvating water, whereas the number of H-bonds between the porphyra-334 and the surrounding water is basically constant all along the energy release process. Occasionally, these H-bonds can fluctuate and their increase, although transitory, favor further the release process. These results provide a clear microscopic picture of the processes responsible for the de-excitation of porphyra-334 and represent a detailed answer to the question concerning the way in which seaweeds get rid of the excess energy without provoking neither their own destabilization not the dissociation of water molecules present in their natural environment.

Figure 1 :

 1 Figure1: The eight configurations of the porphyra-334 molecule used for ground state geometry optimizations and related total energies (a.u.). The label "Unstable" means that proton transfer processes were observed during the geometry optimization procedure. The configuration labeled as 3 can switch to configuration 2, as well as configuration 4 can change to configuration 8 upon optimization. The carbon atoms marked with asterisks inside the cyclohexane rings are linked to three atoms and carry either a double (1.33 Å) or a single (1.54 Å) bond.

Figure 2 :

 2 Figure 2: Schematic picture of the ground (S0) and excited (S1) states of the porphyra-334 molecule. The transition from C to A indicates the initial vertical transition from the ground to the excited state. The transition from A to B indicates the structural relaxation at the excited state. At the point B, the system reverts to the ground state via a nonradiative transition. The transition from B to C indicates the structural relaxation at the ground state. At the point C, the porphyra-334 acquires an excess of kinetic energy as a response to the adsorption of photons of wavelength of 334 nm. Our dynamical simulations describe the process occurring for the transition from B to C and the subsequent relaxation in the range around point C on the ground state.

Figure 3 :

 3 Figure 3: Schematic pictures of the four intermolecular motions.

Figure 4 :

 4 Figure 4: Optimized geometries of configurations 2 and 8. The labeling refers to Figure 1. The total energies of these two structures are separated by only 0.15 eV and they represent the most stable minima. Oxygen, nitrogen, carbon, and hydrogen atoms are represented by red, blue, gray, and white spheres, respectively.

Figure 5 :

 5 Figure 5: Superimposition of the ground (red) and excited (blue) state geometries of the porphyra-334 molecule for the two most stable configurations. The root mean square deviations between the ground and the excited states are 0.51 Å and 0.52 Å for the configurations 2 and 8, respectively.

Figure 6 :

 6 Figure6: Results of the Fourier transform of the velocity autocorrelation functions for the equilibrated part of the simulations. The red line shows the spectrum of the total kinetic energy of the solvent water molecules. The green line show the contribution due to the translational motion. The blue line, instead, is the contributions due to both rotational and vibrational motions with the two peaks at 1600 cm -1 and 3200 cm -1 being the H 2 O molecular bending and stretching, respectively.

Figure 7 :Figure 8 :

 78 Figure 7: Evolution in time of the frequency dependent kinetic energy contribution (a) and analogous time evolution of the difference between non-equilibrated and equilibrated frequency dependent kinetic energy contributions (b) as expressed in equation 2. These results refer to the analysis of the dynamical evolution of the configuration 2. The results of the configuration 8 are basically identical.

Figure 9 :

 9 Figure 9: Evolution of the hydrogen bond network around the porphyra-334 molecule solvate din water for the configuration 2 structure. Panel (a) shows the hydrogen bonding of the solute with solvating water at the beginning of the simulation, whereas panel (b) shows the dynamical evolution at the end of the run where the equilibrium is attained. Red, blue, gray, and black balls and sticks indicate oxygen, nitrogen, carbon, and hydrogen, respectively, while hydrogen bonds are shown as dashed blue lines.

Figure 10 :

 10 Figure 10: Time evolution of the number of hydrogen bonds between water molecules averaged over all the seven independent simulations. The purple lines indicate the evolution in time of the total numbers of hydrogen bonds while the green lines represent the rolling averages. (a) and (b) are results of simulations using configuration 2 and 8, respectively.

Figure 11 :

 11 Figure 11: Time evolution of the numbers of the hydrogen bonds between porphyra-334 and surrounding water molecules averaged over all the seven independent simulations. The purple lines indicate the total instantaneous number of hydrogen bonds, whereas the green lines are the rolling average computed from these time evolutions. (a) and (b) are results of simulations using configuration 2 and 8, respectively.

Acknowledgement

We are deeply indebted to Prof. Shinji Saito at Institute for Molecular Science (Japan) for fruitful discussions and advices. This research was supported by JSPS KAKENHI (No. 25288012), Elements Strategy Initiative to Form Core Research Center (since 2012), and by MEXT, Japan.

We thank computer facilities of ISSP and ITC at The University of Tokyo, RCCS at Okazaki, Japan, and the K computer (Project ID: hp160046) at AICS. A part of the calculations was performed on the HOKUSAI GreatWave system in RIKEN. M.B. is grateful to Pôle HPC and Equipex Equip@Meso at the University of Strasbourg, Grand Equipement National de Calcul Intensif (GENCI) under allocation DARI-A2 A0020906092.

Supporting Information Available

A listing of the contents of each file supplied as Supporting Information should be included. For instructions on what should be included in the Supporting Information as well as how to prepare this material for publications, refer to the journal's Instructions for Authors.

The following files are available free of charge.

• Filename: brief description

This material is available free of charge via the Internet at http://pubs.acs.org/.


Takeshi Baba

Mauro Boero

Katsumasa Kamiya

Hiroyuki Ando

Seiji Negoro

Masayoshi Nakano

Yasuteru Shigeta

Unraveling the degradation of artificial amide bonds in nylon oligomer hydrolase: from induced-fit to acylation processes

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

Introduction

Nylon-6 is a synthetic fiber with a wealth of applications from textiles to electronics industry and, as such, increasingly produced by chemical companies worldwide. The fundamental reaction exploited is the ring-opening polymerization of ε-caprolactams. 1,2 However, this has a non-negligible environmental impact, since it amounts to industrial-scale production of non-biological amide bonds in long fibers containing more than 100 units of 6-aminohexanoate (Ahx). In this respect, an unprecedented discovery dating back to the 1940's is the existence of Nylon-oligomer hydrolases, found in the Arthrobacter sp. KI72, able to degrade the nylon-6. This discovery brought to the attention of the scientific community that this bacterium was able to adapt to the presence of these non-biological compounds and, through natural selection, to develop a new functionality suitable to degrade the artificial Nylon-6 fibers at least for a relatively short period.

Nowadays, three different types of the Nylon-oligomer hydrolases have been identified, namely Ahx-cyclic-dimer hydrolase (NylA), Ahx-liner-dimer (Ald) hydrolase (NylB), and endo-type Ahx-oligomer hydrolase (NylC). [3][4][5][6][7][8] Each one of these hydrolase has its own inherent substrate specificity. Specifically, NylA degrades the amide bond present in a cyclic Ahx dimer to produce a 6-aminohexanoate-liner-dimer (Ald). NylB degrades the amide bond formed in a linear Ahx oligomer, especially Ald, from the N-terminal of the oligomer. NylC degrades amide bonds of both cyclic and linear oligomers. Moreover, NylC is also able to revert nylon-6 into Ahx oligomers. 7,8 Among these three Nylon-oligomer hydrolases, NylB has received special attention in a series of biochemical experiments. According to the reported genetic analyses, NylB has a closely resembling enzyme, NylB * , characterized by a 88% sequence identity. Nonetheless, despite this close similarity, the reactivity to the substrate (Ald) of NylB * is lower than that of NylB by about a factor of 1/200. 3 In an attempt at enhancing the reactivity and to overcome the difficulty in their crystallization, hybrid enzymes of both NylB and NylB * were made by replacing some amino acid residues in NylB * . 9 This resulted in a set of 24 different structures labeled from Hyb-1 to Hyb-24 which were subsequently structurally analyzed by X-ray crystallography or crystallization. A noticeable example of this extensive investigation is represented by Hyb-24; this specific hybrid enzyme contains five amino acid residue substitutions with respect to the pristine NylB * : T3A, P4R, T5S, S8Q, and D15G. Analogously, in the recently obtained Hyb-24DN and Hyb-24DNY a double mutation (G181D and H266N) and a triple mutation (G181D, H266N, and D370Y) with respect to Hyb-24 were realized by one of us. [10][11][12][13] The specific mutational sites were deduced by a comparison of the amino acid sequences of NylB and NylB * (Asp181 and Asn266 in NylB, as opposed to Gly181 and His266 in NylB * ). 9,10 It is also worthy to note that polymerase chain reaction (PRC) applied in random mutagenesis experiments, accidentally resulted in the discovery of the D370Y mutation mentioned above; nonetheless, its role was a matter of debate. Only later 12,13 kinetics analyses allowed to establish that the actual role of D370Y is a stabilization of enzyme substrate complex. At the present state of knowledge, Hyb-24DNY possesses the highest activity among all the known mutants. 12,13 For this reason, hereafter, we shall assume Hyb-24DNY as our reference system and indicate it as the "wild type" (WT) for simplicity. Judging from the similarities in the overall structure and the amino acid sequence, NylB belongs to the class C of the β-lactamase family, which is characterized by the presence of a catalytic triad, Ser, Tyr, and Lys, having the ability to promote the hydrolysis of amide and/or ester bonds. 9 According to X-ray crystallographic analyses, 11,13 the WT can assume two different conformations, (i) a substrate-free form (open form) and (ii) a substrate-bound form (closed from) (see Fig. 1 (a) for details). At variance with the other enzymes of the same family, the distinctive feature of NylB is the presence of a unique amino acid residue, Tyr170, near the catalytic triad, which forms a hydrogen bond (H-bond) with Ald, only in the closed form as sketched in Fig. 1

(b).

Site-directed mutagenesis investigations were conducted on Hyb-24DN, the mutant in which the Tyr170 residue is replaced by a phenylalanine. These studies revealed that the hydrolytic activity of Ald is nearly suppressed, becoming only about 1.4%, and its kinetics parameter values 14,15 , namely the Michaelis constant ( K M ) and the rate constant ( k cat ) could not be determined because of such a low activity. 11 This is an additional difficulty with respect to other mutants such as, for instance, Hyb-24Y (one mutation of D370Y from Hyb-24), Hyb-24DY (double mutation of G181D and D370Y from Hyb-24), Hyb-24DN, and Hyb-24DNY, for which the relevant kinetic parameters could be determined. 13 This can again be ascribed to the fact that the Tyr170 residue is present both in Hyb-24 (NylB * type) and in the WT (NylB type). On these grounds, it can be inferred that Tyr170 is a crucial actor in the hydrolytic activity promoting the cleavage of non-biological amide bonds.

Yet, the atomic-level elucidation of the underlying mechanism and the actual role of Tyr170 present still unclear issues. Part of these unresolved points could be recently addressed by our group 16 by resorting on hybrid quantum mechanics / molecular mechanics (QM/MM) dynamical simulations complemented with free energy sampling methods, allowing to unravel the reaction mechanism in NylB and addressing one of the possible roles of Tyr170 during the acylation process. The latter is the rate limiting step of the series of enzymatic reactions operated by NylB. We also noticed that the residue Asp181 does not contribute directly to the acylation process, because of its relatively distance from the active site (Fig. 1 (c)). Despite being a considerable step forward in the understanding of the uncommon enzymatic activity of NylB, several points remain obscure and a comprehensive picture of the whole mechanism leading to the degradation of amide bonds of non-biological origin has still to be reached. Indeed, several processes are involved, which can be schematically summarized in four principal steps: (1) substrate binding, (2) induced-fit, (3) catalytic (enzymatic) reaction, and (4) final detachment of the substrate. Since the introduction of mutations is expected to affect several processes at the same time, it is necessary to analyze at least the steps (2) and (3), representing the main core of the overall reaction and involving directly the amino acid residues present in the pristine system or substituted upon the mutations. What is left out in the present study are then the initial docking (1) and the final release of the substrate (4), which, however, do not involve chemical reactions. On the other hand, we focus here (i) on the induced-fit of a specific loop region in NylB to the substrate and (ii) on the acylation processes of both the WT and the Y170F mutant. Our purpose is to shed some light into the peculiar hydrolytic activity of NylB and to address the open question about the absence of such an activity in other peptidases and lactamases. Coming to the issue of the docking of the enzyme onto the substrate, we pay special attention here to the mutational effects of D181G on the stabilization of the enzyme-substrate binding aimed at the enhancement of the substrate specificity in comparison with NylB and NylB * . This paper is organized as follows. In section 2 we review the computational approaches used for the reactive QM/MM Car-Parrinello dynamics (CPMD), 17,18 for the Parallel Cascade Selection MD (PaCS-MD) scheme, 19,20 and for the inter-fragment interaction energy (IFIE) analyses based on the fragment molecular orbital (FMO) method. [21][22][23] In the Section 3, we discuss in detail the results of each simulation performed and, on the basis of structural homology analyses, we propose a computationally supported answer to the question concerning the origin of the hydrolytic activity of NylB. Finally, in Sec. 4 we draw our conclusions.

Computational details 2.1 Preparation of target systems

We started the preparation of the prototype structure for the WT enzyme-substrate complex from the X-ray coordinates available at the Protein Data Bank (PDB entry: 2zma). In this system the Ser112 residue is replaced by Ala to crystallize the complex. To revert back to the system targeted in this study, we replaced Ala112 with the original Ser112 residue and added all the missing amino acid residues to obtain an appropriate WT model. Analogously, the structure of the Y170F mutant was reconstructed by replacing Tyr170 with Phe170 in the WT modeled as explained above, whereas the Hyb-24 was obtained by substituting Asp181, Asn266, and Tyr370 with Gly181, His266, and Asp370 in WT, respectively. The target WT system used in the simulations consists of the above enzyme-substrate complex solvated in water (described by a TIP3P force field) 24 plus Na + counter-ions to neutralize the system; such a model amounts to 41403 atoms.

The system was prepared by first equilibrating the solvated complex via classical molecular dynamics (MD) for about 10ns in an NPT ensemble, by setting the target temperature and pressure at 300 K and 1 atm, respectively. To this aim, the PMEMD module (SPFP precision model) of the AMBER 25,26 program suite was used. The force field adopted for the amino acid residues was a standard Amber ff99, while for Ald the force field was generated ad hoc with the AMBER GAFF tool (see Ref. 27) for details using the Antechamber program of the AMBER package. 25 The solvated system is placed in a simulation cell of 77.57 x 76.62 x 69.30 Å 3 where periodic boundary condition are applied and the SHAKE algorithm 28 was used all across MD simulations to satisfy bond geometry constraints. The long-range interactions are evaluated by means of the particle mesh Ewald method with a cutoff of 12.0 Å. 29 The protonation states adopted here for each amino acid residue are the ones originally determined by the group of J. H. Jensen [30][31][32][33] and partly the ones deduced by direct inspection of their local hydrogen bonding environments. In particular, Lys115, which is one of the fundamental amino acid residues in the active site, was set to be deprotonated (-NH2) to promote the hydrolysis.

Reactive QM/MM Car-Parrinello simulations for the acylation

The acylation reaction, which involves cleavage and formation of chemical bonds, was simulating within a reactive hybrid QM/MM approach 34,35 complemented with metadynamics [36][37][38][39][40] in which the QM driver is the Car-Parrinello molecular dynamics (CPMD) 41 based on the density functional theory. 42,43 The aim is the investigation of the reaction mechanism of the Y170F mutant and the estimation of the free energy difference between WT and Y170F. This approach has already been extensively used and benchmarked many biological systems. [44][45][46][47][48][49][50][51] In our specific systems, the embedded QM region is composed of 103 atoms in the case of the WT and 102 atoms in the case of the Y170F mutant (Fig. 1 b). The total net charge for the QM subsystem equal to -1 e and monovalent hydrogen-like link atoms are used to saturate dangling bonds at the interface between the QM and the MM regions. Further details concerning the QM subsystems can be found in the supporting information (SI-1). The CPMD package 52 was used for all these simulations. The exchange correlation functional used here is the one proposed by Hamprecht, Cohen Tozer and Handy (HCTH), 53 the core-valence interaction is described by Martins-Troullier pseudopotentials 54 and valence electrons are represented in a plane wave basis set with an energy cut off of 70 Ry. All QM/MM CPMD simulations were done in a canonical NVT ensemble at a target temperature of 300 K. The fictitious electron mass for the integration of the CPMD electronic equations of motion was set to 400 a.u. and the time step adopted was 4 a.u. (0.096 fs). Classical MM interactions were described by an Amber ff99 force field. Before moving to the reactive metadynamics 36,37 simulations, the system was equilibrated via classical MD in an NVT ensemble set-up for 1 ns, then, subsequently, and unbiased QM/MM CPMD dynamics of the same type was performed to equilibrate also the newly added electronic degrees of freedom for about 2 ps. After the equilibration, metadynamics simulations were done with the purpose of understanding whether or not the role of the catalytic triad Ser112-Lys115-Tyr215 is identical both in the WT and in the mutant. Although the reaction mechanism for the WT has already been reported elsewhere, 16 we repeated this simulation using one single collective variable as explained later. On one hand, this provides a further and independent check of our former results and, on the other hand, the identical scheme applied to the Y170F mutant allows for a direct comparison of the relative reactivity of the two systems. To the benefit of the explanation of the parameters adopted, we recall that, in its Lagrangean formulation, the metadynamics scheme is summarized by the extension of the Car-Parrinello Lagrangean L CP to include the selected reaction coordinates, referred to as collective variables (CVs) hereafter, namely ,

Concerning the single CV, s, used in the present case, we selected the distance between O g atom (OG) of Ser112 and carbonyl carbon (C1) of Ald393 (indicated by a green arrow in Fig. 1 c), which properly represents the nucleophilic attack of Ser112 to the amide carbon present in the substrate.

The history dependent potential V t, s a ( ) was constructed by accumulation of small Gaussian penalty functions having height width of 0.5 kcal/mol and 0.3 au and a new Gaussian function was added every 140 steps (13.44 fs), which will be referred to as one meta-step hereafter (Details of other parameters are shown in SI-2). The history-dependent term in principle converges asymptotically to the free energy surface (FES) F s a ( ) which can then be computed as the average

( )     simul diff , 1 ) ( diff simul t t dt t s V t t s F a a , (3)
to reduce statistical errors 55 , tsimul being the total simulation time and tdiff the time at which CVs start diffusing in all the sampled phase space. For the present purpose, this can be assumed to be the completion of the acylation. Indeed, in this case, we followed the simulation until the first reaction occurs, since we focus specifically on the relative free energy barriers in WT and Y170F for this specific process. Because of the specific choice of collective variable, the system is unable to realize the reverse reaction unless reflecting walls able to reverse the momenta of the leaving atoms are used (as e.g. in Ref. 45). Since this would result in a lack of momentum conservation, we avoided it to prevent any further biasing. However, this procedure is sufficient to estimate, at least approximately, free energy differences. Results were analyzed with our own program and visualized by VMD. 56

PaCS-MD simulation for induced-fit process

To investigate the difference between the WT and the Y170F mutant in the induced-fit processes, we resorted on the PaCS-MD approach, which generates reactive trajectories from a reactant (open form) to a product (closed form) by combining a series of short-time simulations. 19,20 Each one of the two structures was extracted from a 10 ns classical MD simulation. The induced-fit mechanism concerning the WT has already been discussed in a former work 20 and we refer to those results for the ongoing discussion. Instead, for the Y170F mutant PaCS-MD was done for 25 cycles starting from 10 different initial configurations sampled on the basis of the same selection rule formerly adopted, specifically the root mean square displacement (RMSD) of the loop-segment (Asn166-Val177) except for all hydrogen atoms, invisible to X-ray crystallographic inspection, and Cα atoms of the 111 th to 116 th residues. This specific selection rule has been shown to efficiently trigger the induced-fit for WT. 20 To evaluate the free energy change due to the induced-fit, we adopted the weighted histogram analysis method (WHAM) [START_REF] Hankins | [END_REF] , where 2000 reference structures (200 frames randomly selected selection in each one of 10 independent reactive trajectories) were used for the WHAM inspection. A biased umbrella potential with a coupling parameter k = 0.5 kcal/mol Å 2 was imposed on the loop-segment, apart from H atoms, the nitrogen atom of the amide bond of Ald, and the Cα atoms of the 111 th to 116 th residues. The free energy obtained in this way was then projected onto two specific distances: The separation between C d atom (CD) of Glu168 and Cα atom (CA) of Ser217, d1, and the one between H z atom (HH) of Phe170 (replacing O z atom (OH) of Tyr170 in the WT) and nitrogen atom in amide bond (N3) of Ald393, d2, which clearly correspond to the reaction coordinates previously adopted for the WT. 20

FMO calculations of enzyme substrate bound state of Hyb-24 and Y170F

To better quantify the interaction between ALD and the specific amino acid residues involved, namely Asp181 in Hyb-24 and Tyr170 in the Y170F mutant, we made use of an inter-fragment interaction energy (IFIE) analyses based on the fragment molecular orbital (FMO) method. 21 The total energy, E, in FMO framework is expressed as

E  E I I å  E IJ  E I  E J ( ) I >J å , (4)
where EI and EIJ represent the fragment energies of I-th monomer and IJ-th dimer, respectively. 53 The IFIE value is defined as 59

DE IJ  E IJ '  E I '  E J ' + Tr ∆ P IJ V IJ [ ] , (5)
where ∆ P IJ is the matrix constructed by computing the difference between the IJ-th dimer and two I-th and J-th monomers, while V IJ is the environmental electrostatic potential of the IJ-th dimer with respect to all the other remaining fragments. The superscript on E ' indicates the bare fragment energy without the inclusion of the environmental electrostatic potential.

All FMO calculations were carried out as implemented in the ABINIT-MP program suite 60 at the FMO2 MP2/6-31G(d) level. To reduce the computational costs, we limited the calculations to one single snapshot extracted from a 10 ns trajectory obtained upon classical MD. Moreover, explicit solvent effect on the stability of the enzyme-substrate bound form was accounted for by the inclusion of all water molecules included in a region of 15 Å around the substrate. 61,62 The results obtained by IFIE were visualized with the BioStation Viewer. 63 Tyr215 protonates the nitrogen atom of the amide bond (r3~1.0 Å). Such a reaction mechanism for the Y170F mutant is consistent with that of WT, 16 (see SI-3) and the estimated free energy barrier is reported in Table 1. We recall that a typical error bar of 1-2 kcal/mol affects these values. Errors can be kept under control in the variational formulation of metadynamics used in our approach, which allows collective variables to oscillate around the "exact" solution in a way identical to the one in which the Car-Parrinello fictitious electron kinetic energy fluctuates around the exact BO trajectory.

Results &Discussion

The dynamics of these extra degrees of freedom can be explicitly controlled by tuning the fictitious mass (kinetic term) and spring constant (harmonic term) and/or by suitable thermostats 55 and, thus, can be smoothened. Nothing, instead, can be done for the underlying DFT level. 83 The most remarkable feature is the larger value of this barrier for the Y170F mutant with respect to the WT, namely 11 kcal/mol above the analogous barrier for the WT, 16 in qualitatively good agreement with experiments in which mutations of Y170F resulted in a dramatic decrease in the Ald degradation activity to about 1.4% with respect to WT. 11 Table 1 Fig. 3 To elucidate the differences between the WT and its Y170F mutant during the acylation, the structural changes, especially when Tyr215 provides the H atom to the amide N atom of the substrate, we compared the evolution of the mutant with respect to the WT system as sketched in the panels of Fig. 3. Before the acylation, the amide bond and the neighboring atoms are all located on an amide plane. Then, when the acylation occurs, the O atom belonging to the side chain of Ser112 attacks the amide bond forming an sp 3 bonding configuration which includes the amide carbonyl C atom, resulting in the breaking of the initial planar symmetry mentioned above.

Simultaneously, the amide N atom also assumes an sp 3 . In addition, it is expected that hydrogen bond might contribute the acceptance tendency of hydrogen atom from Tyr215 as well as the control of direction of NH group considering the Fig. 3 (a). This indicates rather unambiguously that the formation of an H-bond to the amide NH group is indispensable to trigger the H transfer from Tyr215 and that Tyr170 provides the fine control needed to keep the local geometry in a condition suitable to give rise to an efficient breaking of the amide bond. We can then infer that Tyr170 is the fundamental additional residue that must complement the catalytic triad to realize the acylation in this enzyme.

Its absence, as in the case of the mutant, results clearly in the suppression of this specific activity or in a larger activation barrier that makes acylation highly improbable. For amines, two different sp 3 configurations are in general allowed, as shown in Fig. 3-(e), and the conformational change between them is well known as an umbrella inversion, whose barrier is about 6 kcal/mol for an ammonia molecule in the gas phase. 65 By reasoning on these grounds, the order of magnitude of the free energy difference between WT and the Y170F mutant in terms of activation barrier (~11 kcal/mol, according to our estimates) might be compatible with the results obtained via the umbrella inversion technique also in a protein environment.

Comparing our result with theoretical studies on general serine proteases, [66][67][68] structures like the one reported in Fig. 3-(e) were found, for instance, in Fig. 4 of Ref. 62 and discussed in the context of the stereochemistry of tetrahedral intermediates in Ref. 61. In the typical proteases, such as, for instance, chymotrypsin and trypsin, a substrate pocket able to bind to their specific substrate is always present. In the substrate pocket, the chymotrypsin and the trypsin recognize bulky aromatics and side chain in a peptide, respectively, enabling them to identify the site to be cleaved. Instead, the substrate-binding site of NylB is more like a cleft, and in principle can bind substrates so weakly that simple mechanical fluctuations in substrate and changes in the organization of water molecules around the substrates result in a decrease of the reactivity. Only after the induced-fit of Tyr170, allowing to smooth the fluctuations of the substrate in the vicinity of the catalytic triad and to prevent approaching of water molecules during the acylation phase, NylB can degrade the Nylon-oligomers.

So far, focusing on the acylation processes, we established the different behavior of the WT and the Y170F mutant, providing a solid argument for the decrease in the hydrolytic activity of the latter.

Nonetheless, a still unaddressed question is the actual role of the Tyr170 residue in the induced-fit to bind the substrate.

(A) and (B), but also the other remaining two, are rather wide and characterized with not very steep slopes with respect to the WT. Although Glu168 is a residue common to both the WT and the mutant, in the first case it forms an H-bond with the substrate and Ser217; conversely fluctuation of Phe170 in the Y170F mutant prevents Glu168 from forming any H-bond with Ser217. Despite this important structural difference, the highest free energy difference in the FES of the Y170F mutant is nearly identical to that of the WT. These findings indicate that the mutational (structural) effect due to the induced-fit on the FES landscape is remarkable, while the stabilization energy is small from a kinetic point of view.

Substrate binding stability of Hyb-24 and Y170F mutants

In an attempt at quantifying the effect of each residue on the overall stability of the enzyme-substrate bound form, such as Hyb-24, in particular the so-called D181G mutant, and the Y170F mutant, we performed IFIE analyses based on the FMO method. More precisely, we focused on the IFIE values between each ALD fragment and the most relevant amino acid residues, such as Tyr170 and Asp181 (see Fig. 5). In order to analyze contributions of IFIE in more detail, we further divided the substrate into three distinct fragments: the C-terminal (Ald(1)), the amide group (Ald(2)), and the N-terminal(Ald(3)), as shown in Fig 5 . The result of this analysis is summarized in Table2.

The IFIE values for the WT have been already elsewhere 56 and shown in SI-5. According to Table 2, Asp 181 is characterized by an attractive interaction with Ald(3) stronger than that of Tyr170 with Ald(2). Conversely, in Hyb-24 the interaction between NylB and ALD(3) is lacking. Although the IFIE value to Ald(2) of Phe170 in the Y170F mutant becomes weaker than that of Tyr170 in the WT, the decrease of the IFIE is less pronounced than that to Ald(3) of the Hyb-24. These results suggest that Asp181 is quantitatively, energetically speaking, an amino acid residue more important than Tyr 170 for the specific substrate binding. Indeed, available kinetics data on Hyb-24DY, resulting from a double mutation of G181D and D370Y from Hyb-24,

( k cat K M  0.35 s 1 mM 1 , K M  7.1 ±1.28mM 1
) and of Hyb-24Y, in which a single mutation of D370Y from Hyb-24 is present, ( k cat K M  0.016 s 1 nM 1 , K M  39.1 ± 6.73mM 1 ), 13 indicate rather clearly that Asp181 contributes to enhance the substrate specificity.

Homology analyses on NylB and possible molecular evolution route for Ald hydrolytic activity

Our computational results, corroborated by experiments, suggest that, from the evolutional point of view, the introduction of Y170 in the loop region is at the origin of the acquisition of the hydrolytic activity from NylB. To bring support to our assumption, we completed our study by investigating any homologous enzyme exists possessing a similar evolution for the substrate binding. To this purpose, we performed structural homology searches for the WT (PDB id: 2zma) by using the DALI program 69 and selected structures preserving the local configuration around Tyr170 similar to NylB. This search, performed by using MultiSeq2.0 with the STAMP algorithm. [70][71][72] , brought to our attention the EstU1 system, having a high homology to NylB, as shown in Fig. 6 (a). The specific substrate of EstU1 is a cephalothin, which has bulky functional groups. Around the binding site of the substrate, EstU1

presents three amino acid residues, Tyr172, Phe174, and Phe190, which correspond to Glu168, Tyr170, and Asp181 in NylB, respectively. These three aromatic amino acid residues recognize the thiophene ring of cephalothin cooperatively by hydrophobic interactions (see Fig. 6 (b)). The remarkable difference between NylB and EstU1 is the location of the helix structure and the loop-segment evidenced in Fig. 6. In addition to the similarity in the structure, the side chain of Phe174 in the loop-segment of EstU1 undergoes a conformational change during the substrate binding process. This property is resembles closely the one of NylB. We can then infer that, upon acquiring Tyr170, NylB might have undergone a natural selection process to adapt to a nylon-oligomer-rich environment by skillfully using the loop region for the induced-fit to the substrate and by catalyzing the acylation via a decrease if the activation barrier.

Fig. 6 Towards a rational design of enzymatic activities

Having unraveled the role of Tyr170 in the WT, which contributes mainly to the acylation step to regulate the direction of a hydrogen bond, we have now a convincing set of results allowing to proceed further. All our results concur to support the notion that NylB * needs not only the catalytic triad Ser-Tyr-Lys, like the class C β-lactamase family, but also Tyr170 to degrade non-biological amide bonds present in nylon-oligomers. In particular, acquiring Tyr170 in NylB * is the key factor for triggering the reaction by decreasing the activation barrier, in other words by decreasing the rate constant k cat . This is confirmed by a comparison with Hyb-24 (nearly identical to NylB * ), which possesses as well a weak reactivity to Ald. In the change from NylB * to NylB, Asp181 provides the required substrate specificity, resulting in the decrease of k cat .

Recent advances in computer facilities and related theoretical methods allow nowadays for a de novo design of functional proteins. 73 In the last decade, researchers working in this field have succeeded in producing new enzymes called "theozymes", which were entirely designed on a computer by inserting a suitably selected transition state structure into protein scaffolds. [74][START_REF] Privett | Proc. Natl. Acad. Sci[END_REF][START_REF] Mills | [END_REF][77] Unfortunately, for these theozymes, the resulting catalytic proficiency, defined as

Catalytic proficiency  k cat / K M ( ) k uncat in M 1 ( ) , (6)
was still inferior to the natural one. In this expression k uncat is the rate constant in the absence of the enzyme, 78 while k cat is the rate constant of the theozyme, lower than the natural one. 79 The lesson to be retained is that looking just at the structure of the building block (transition state in this case) to be inserted inside a given scaffold is insufficient to obtain an enhancement of the catalytic activity. In this respect, in our study we have shown that Tyr170 and Asp 181 play an important role in direct and indirect decrease of k cat by stabilizing both the transition state and the substrate binding structure, thus allowing for an efficient degradation of the non-biological amide bond. Compared with the strategy followed for the theozyme, the change realized in the passage from Hyb-24 (NylB * type) to WT or Y170F (NylB type) could resemble the processes of the transition state design and the active site re-design briefly summarized above. Nonetheless, the theozyme was designed and realized in a laboratory and never underwent any natural selection process. Its absence in nature might be ascribed to its relative rigidity, which, contrary to other enzymes, implies a lack of dynamic motions of the flexible loop regions. We can then expect that taking into account also the response of the system to natural selection, as in the case of NylB, can be the key for a rational design of enzymatic reactions with appreciable catalytic proficiency, although we have to acknowledge that a major difficulty is still represented by the prediction of the possible conformational changes an enzyme could assume. 80 Some hints could be obtained by the development of integrated schemes making use of coarse grained MD 81 or Fluctuation Flooding Method 82 jointly with quantum chemical transition state design not too far from the strategy adopted for the theozyme.

Conclusions

A combination of different theoretical methods was used to inspect in detail a series of enzymatic reaction representing the fundamental stages in the degradation of non-biological amide bonds operated by the nylon oligomer hydrolase NylB. To this aim a comparative study of the WT enzyme and its closely related mutant was conducted, unveiling the crucial role of two evolutionally important amino acid residues, Tyr170 and Asp181, in (1) the enzymatic reactions, (2) the induced-fit process, and (3) the stability of the enzyme-substrate binding. We found at which steps of this complex process Tyr170 and Asp181 are involved and which are their respective net effects on the overall process. Specifically, Tyr170 mainly contributes to the acylation step and plays a subsidiary role in the stabilization of the enzyme-substrate complex. According to the free energy changes during the acylation and the induced-fit steps upon the mutation, in the mutant the former process requires the overcoming of an activation barrier which is 11 kcal/mol higher with respect to the WT case. This is clearly a stumbling block to acylation that makes the reaction very unlikely for the Y170F system. Conversely, for the induced-fit process of the Y170F mutant, although the free energy landscape changes drastically with respect to the WT case, the free energy barrier is nearly unchanged, suggesting that the kinetics of the induced-fit is barely affected by the mutation or, at least, this mutation. Although the inter-fragment interaction energy (IFIE) of Phe170 in the Y170F mutant is lower than that of Try170 in the WT, such a decrement is smaller than that of Hyb-24.

Based on all these results, we can conclude that Tyr170 is the essential crucial amino acid residue to promote the degradation of non-biological amide bonds by forming a novel type of a catalytic quad, as opposed to standard catalytic triads, responsible for both the stabilization of the transition state and the restraining of fluctuations in the substrate binding structure. On the other hand, Asp181 mainly contributes to the stabilization of the enzyme-substrate complex (substrate specificity), while the influence on the hydrolytic activity is negligible. Each substrate is shown as yellow sticks.

Table 1. Maximum free energy differences in the acylation step for the WT and the Y170F mutant.

We recall that typical error bars of 1-2 kcal/mol from the underlying level of DFT adopted affect

Table of contents

To unravel what provides the ability to degrade non-biological amide bond with Nylon oligomer hydrolase, we investigate the process from induced-fit to acylation by a combination of different theoretical methods.

Fig. 1

 1 Fig.1

Fig. 2

 2 Fig.2

  geometry. During this process, in the case of the WT, the Tyr170 contributes to the reaction with its own -OH group: This -OH forms one H-bond with the -NH moiety of the amide bond of Ald and this, in turn, forces the nitrogen atom in Ald to keep its lone pair oriented toward Tyr215. In such a configuration, the amide N pointing toward Tyr215 is ready to accept the H atom from Tyr215, as confirmed by the formation of a new and tight H-bond between the phenol OH and the amide NH group. At variance with the WT case, in the Y170F mutant (Fig.3-(b)) the absence of Tyr170 and its H-bond does not allow the amide NH group to assume a configuration suitable to accept the H atom from Tyr215, thus suppressing the amide bond cleavage activity. This picture is perfectly consistent with the fact, recently demonstrated64 , that the specific orientation of the lone pair of the nitrogen of the scissile amide bond during catalysis is a crucial feature of enzyme catalyzed amide bond hydrolysis. After a few additional meta-steps, a water molecule from the surrounding solvent approaches to form an H-bond with the amide NH group as shown in Fig3-(c), thus replacing the role of Tyr170 in the WT. From the point of view of the local H-bond environment only, this conformation resembles closely the one of the WT shown (Fig3-(d))

Fig. 5

 5 Fig.5

Figure 1 .

 1 Figure 1. Structure and main amino acid residues of NylB. (a) Induced-fit process: Superposition of the NylB structures in the open form (black) and in the closed form (yellow). A large movement of the loop-segment N166-V177 toward the substrate occurs, indicated by the red arrow enclosed in the dotted red circle. (b) The QM region, located inside the red circle in (a), used in our hybrid QM/MM simulations. QM atoms are shown as balls and sticks, whose colors code, here and in the all the following figures, are gray for C, red for O, blue for N and white for H, and the yellow ribbon shows the secondary structure. (c) The reaction coordinate (red arrow) and other fundamental distances (green, blue and magenta) used to drive and monitor the acylation reaction. Details are given in the text and in Fig. 2.

Figure 2 .

 2 Figure 2. Evolution of the reaction coordinate r1 and the most relevant distances during the metadynamics simulation for the Y170F mutant: r1 = N3(Ald393)-C1(Ald393) (red), r2 = OG(Ser112)-C1(Ald393) (green), r3 = HH(Tyr215)-N3(Ald393) (blue), and r4 = HG(Ser112)-NZ(Lys115) (magenta). Panel (a) shows whole simulation, while panels (b), (c), and (d) are the magnifications of the three crucial events inside the dashed circles in (a).

Figure 3 .

 3 Figure 3. Local structures and related differences between Y170F and WT relevant to the acylation process. Panels (a), (b), and (c) correspond, respectively, to the events (b), (c), and (d) in Fig. 2. Panel (d) shows the WT for comparison. A schematic picture of the difference between Y170F and WT is shown in (e).

Figure 4 .

 4 Figure 4. Free energy landscapes for the induced-fit process in the Y170F mutant(a) and the WT (b), according to the results of Ref. 20.

Figure 5 .

 5 Figure 5. (a) Schematic view of the interactions between the most relevant amino acid residues and the substrate. For clarity, the ALD is divided into three fragments, the C-terminal, the peptide site, and the N-terminal, labeled as ALD(1) (blue), ALD(2) (pink), and ALD(3) (green), respectively. These fragments are the ones used in the FMO analysis shown in panel (b).

Figure 6 .

 6 Figure 6. Result of the structural homology search (a) and details of the surrounding substrate (b) Using the DALI program and selecting only structures preserving the local environment around Tyr170 in WT (pdb id 2ZMA), EstU1 (pdb id 4IVK) was discovered. For the structural aliment, we adopted the STAMP algorithm. The corresponding analogous main regions between WT and EstU1 are shown in the two adjacent panels in (a). The semi-transparent elliptical region on the top left of each structure indicates the two segments of 4IVK, and 175-185 of 2ZMA, respectively.

  kcal / mol

Figure 2 TFigure 3 T

 23 Figure 2 T. Baba et al.

Figure 6 T

 6 Figure 6 T. Baba et al.

Table 2 .

 2 The computed IFIE values between the 170 th and the 181 st amino acid residues and each fragment in Ald.

				IFIE a
	System	Residue
			ALD(1)	ALD(2)	ALD(3)	Total
	WT b	Tyr170	-2.07	-6.01	0.73	-7.35
		Asp181	18.31	-2.63	-115.27	-99.59
		Phe170	-4.69	0.33	-1.44	-5.80
	Y170F
		Asp181	20.19	-3.05	-118.21	-101.07
		Tyr170	-4.87	-4.79	-4.34	-14.00
	Hyb-24
		Gly181	-0.52	0.13	1.22	0.83

a unit of kcal / mol, b From Ref.

59

ACKNOWLEDGMENT We thank Mr. Masanari Matsuoka for insightful discussions in the search for structural homology. This work is supported by Grant-in-Aid for Scientific Research of the Innovative Areas "Computics", "Dynamic ordering and biofunction", and Photosynergetics" (No. 25104716, 26102525, and 26107004), and a Grant-in-Aid for Young Scientists (B) (Grant No. 22740259) from Japan Society for the Promotion of Science (JSPS), and by the Strategic Programs for Innovative Research (SPIRE), MEXT, and the Computational Materials Science Initiative (CMSI), Japan. T.B. expresses his special thanks to JSPS Research Fellowship for Young Scientists (No. 02503029). Theoretical calculations were partly performed at the Research Center for Computational Science, Institute for Molecular Science (Japan) and Center of Computational Materials Science, Institute for Solid State Physics, The University of Tokyo, (Japan). M.B. and T.B. acknowledge the Direction Informatique (Pôle HPC) of the University of Strasbourg for providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project. K.K. acknowledges the Cooperative Research Program of Institute for Protein Research, Osaka University

Induced-fit process

To complete the picture about the overall role of Tyr170, we investigated the induce-fit process.

To this aim, the free energy surface (FES), projected onto two reaction coordinates d1 and d2, was estimated by WHAM making use of a set of reactive trajectories generated via PaCS-MD simulations (Fig. 4 The highest free energy barrier is the one separating (C) from (D) and its value is about 2.4 kcal/mol. Additional information can be found in the supporting information (SI-4). Fig. 4 To elucidate the mutational effects of Y170F on the induced-fit process (Fig. 4 (a)), we compared the present results with the previous ones for the WT. 20 Focusing on the FES landscape of both systems, two major features are worthy of note. The first one is the number of (meta-) stable structures realized by the two systems, namely six in the case of the WT, as opposed to four for the Y170F mutant. The second one is the size of each basin where these local minima are located. In the case of the Y170F mutant, the areas of the free energy landscape around the minima, especially


A Scholtissek
email: scholtissek@stfs.tu-darmstadt.de

P Domingo

L Vervisch

C Hasse

A self-contained progress variable space solution method for thermochemical variables and flame speed in freely-propagating premixed flamelets

Keywords: premixed flamelet, progress variable space, gradient equation, self-contained, laminar flame speed

Flamelet models for premixed combustion, which are based on equations formulated and solved in progress variable space, have been proposed in the past, but have not been adopted for chemistry reduction methods. This is due to one limitation of these models: they need a closure for both the magnitude and the shape of the gradient (or scalar dissipation rate) of the progress variable, which is essential for an accurate prediction of the flame displacement speed. So far, solution methods for the aforementioned models require gradient information as an input, which is either modelled and non-generic, or extracted from a previous physical space flame solution for the analogous problem. The objective of this work is to provide a self-contained solution method for freely-propagating premixed flamelets in progress variable space, by solving an additional flamelet equation for the gradient of the progress variable. With this, the novel method provides both magnitude and shape of the gradient. Studying hydrogen-air and methane-air configurations, it is demonstrated that an accurate prediction of the laminar flame speed without the necessity for further input parameters can be obtained.

Introduction

Flamelet [1] and flamelet-based models, such as Flamelet-Generated Manifolds (FGM) [2], Flame Prolongation of ILDM (FPI) [START_REF] Gicquel | Proc. Combust. Inst[END_REF] or the Reaction/Diffusion Manifold (REDIM) concept [START_REF] Bykov | [END_REF], have been successfully used as reduction methods for CFD-simulations of chemically reacting flows. Thermochemical states of generic flame configurations are computed a-priori and stored in flamelet tables as a function of few control variables, such as mixture fraction and/or progress variable. In comparison to full chemistry solutions, these methods exhibit a dramatically reduced computational cost by solving only for the flow and the controlling variables, while the remaining thermochemical state is obtained from the flamelet table.

Among flamelet models, where equations are formulated and solved in either mixture fraction space (non-premixed combustion) or progress variable space (premixed combustion), only the former have been extensively used and extended for a wide range of conditions and applications. Few works concentrate on premixed flamelet models in composition space [START_REF] Nguyen | Combust[END_REF][6][7], although these models show great potential for the aforementioned reduction methods. This is supported by the many successful applications with their physical space counterpart, the FGM method [2], where one-dimensional flamelet equations are formulated and solved in a flame-adapted coordinate system. There is one major reason why progress variable space flamelet models for premixed combustion have not been widely used: they need a closure for both the magnitude and the shape of the gradient (or scalar dissipation rate) of the progress variable, which is essential for an accurate prediction of the flame displacement speed. In [START_REF] Nguyen | Combust[END_REF] and [6] premixed flamelet equations have been solved, but profiles for the scalar gradient were input variables and had to be modelled or obtained from physical space solutions. In [7] premixed flamelet equations were analyzed and their validity was also demonstrated for turbulent premixed flames based on direct numerical simulations (DNS). Furthermore, the gradient (or scalar dissipation rate) of the progress variable is a central parameter in most modeling approaches for turbulent premixed flames [START_REF] Kolla | Scalar dissipation rate based flamelet modelling of turbulent premixed flames[END_REF] since it is linked to the fundamental quantities in turbulent combustion [START_REF] Veynante | [END_REF].

It is the objective of this work to show how a closure for the scalar gradient of the progress variable can be obtained for unstrained premixed flamelets. Therefore, an additional flamelet equation is solved for the gradient, resulting in a self-contained model which also allows the laminar flame speed to be computed in the progress variable space. Being the counterpart to the freely-propagating flame in physical space, the model presented here covers the most important canonical flame configuration for premixed combustion.

Self-contained Flamelet Formulation

A freely-propagating premixed flame can be described by the following equation set [START_REF] Kee | PREMIX: a fortran program for modeling steady laminar onedimensional premixed flames[END_REF]

∂ρu ∂x = 0 , (1)
ρc p u ∂T ∂x = ∂ ∂x λ ∂T ∂x -ρ N s k c p,k Y k V k ∂T ∂x + ωT , ( 2
)
ρu ∂Y i ∂x = - ∂ ∂x (ρY i V i ) + ωi , (3)
where x is the spatial coordinate of the one-dimensional physical space, ρ is the density, u the flow velocity, c p the specific heat capacity, T the temperature and λ the thermal conductivity. The quantities Y i and V i denote the mass fraction and the diffusion velocity of species i, respectively, and the chemical source terms for temperature and species are represented by ωT and ωi . The mass flux

ρu = ρ 0 s l , (4)
is an eigenvalue for the above equation system and is directly linked to the laminar flame speed s l (ρ 0 denotes the density of the unburnt gases). The progress variable Y c is introduced as Y c = i α i Y i , where α i is a species-specific weighting factor. An analogous summation of the species equations (Eq. ( 3)) yields a balance equation for the progress variable:

ρu ∂Y c ∂x = - ∂ ∂x (ρY c V c ) + ωc , (5)
where the diffusive flux of the progress variable is defined as ρY c V c = ρ i α i Y i V i . Introducing the gradient of the progress variable g c = ∂Y c /∂x and using Eqs. (4) and 5, an expression for the laminar flame speed is obtained

s l = 1 ρ 0 g c - ∂ ∂x (ρY c V c ) + ωc . (6)

Transformation Rules

In the following sections, Eqs. (2), (3) and ( 5) are transformed and the progress variable Y c is chosen as the control variable. The transformation rule from the physical space coordinate system x to the flame-adapted progress variable space Y c reads

∂T ∂x = ∂Y c ∂x ∂T ∂Y c , (7)
which is applied analogously to other scalar quantities.

Flamelet Equations for Equal Diffusivities

Assuming equal diffusivities for species and temperature, or equivalently unity Lewis numbers, the species diffusion velocity is expressed as

V i = - D Y i ∂Y i ∂x , (8)
where the diffusivity is defined as D = λ/(ρc p ). Applying the transformation rules in Eq. ( 7) to Eqs. ( 2) and

(3) (detailed derivations are found in the supplementary material) yields the premixed flamelet equations for temperature and species

ρDg 2 c ∂ 2 T ∂Y c 2 -ωc ∂T ∂Y c + ωT c p + ρDg 2 c c p ∂c p ∂Y c ∂T ∂Y c + ρDg 2 c k c p,k c p ∂Y k ∂Y c ∂T ∂Y c = 0 , (9)
ρDg 2 c ∂ 2 Y i ∂Y c 2 -ωc ∂Y i ∂Y c + ωi = 0 . (10)
As outlined earlier, an additional flamelet equation is solved for g c , which is obtained by differentiating Eq. ( 5)

∂ ∂x u ∂Y c ∂x = - 1 ρ ∂ ∂x (ρY c V c ) + ωc ρ . (11)
Applying the transformation rules to this expression yields

ag c + ωc ρ ∂g c ∂Y c = g 2 c ρ ∂ 2 ∂Y c 2 (ρDg c ) - g 2 c ρ 2 ∂ρ ∂Y c ∂ ∂Y c (ρDg c ) + g c ∂ ∂Y c ω ρ . ( 12
)
where a is the strain rate, defined as [START_REF] Hasse | [END_REF][START_REF] Peters | 5 th US Combustion Meeting[END_REF]

a = - ∂u ∂x . (13)
Using Eq. ( 1), a can be related to the eigenvalue of the flame

∂ρu ∂x = 0 = -ρa + ρ 0 s l ρ ∂ρ ∂x . ( 14
)
Applying the transformation rule in Eq. ( 7) and using Eq. ( 6), the strain rate equals

a = - g c ρ 0 s l ρ 2 ∂ρ ∂Y c = - 1 ρ 2 ∂ρ ∂Y c - ∂ ∂x (ρY c V c ) + ωc . (15)
This shows that a is non-zero also for an unstrained premixed flame, which is due to thermal expansion since it is a function of the density gradient. Introducing Eq. ( 15) into Eq. ( 12) yields

g 2 c ∂ 2 ∂Y c 2 (ρDg c ) -ωc ∂g c ∂Y c + g c ∂ ωc ∂Y c = 0 . ( 16
)
Considering their origin in the transformation, the terms are interpreted as (left to right): diffusion, Y cdrift (transformed convective term) and chemical source term. Equations ( 9), ( 10) and ( 16) represent a closed equation system which can be solved in progress variable space.

Flamelet Equations with Differential Diffusion

In this section, the premixed flamelet equations are generalized to capture also differential diffusion effects. Keeping in mind that the species diffusion velocity V i generally contains gradient information, it can be written as

ρY i V i = ρY i Ṽi ∂Y c ∂x , (17)
where gradients have been transformed according to Eq. ( 7) and Ṽi represents a diffusion velocity in progress variable space. Instead of gradients with respect to physical space coordinates, Ṽi contains gradients with respect to the conditioning variable Y c and can therefore be computed from standard transport libraries on a Y cgrid. Besides unity Lewis number diffusion, a mixtureaveraged diffusion model [START_REF] Curtiss | [END_REF] with a correction velocity [14] for mass conservation is used in this work. In principle the approach allows arbitrary diffusion models to be transferred to progress variable space. The diffusive flux of the progress variable is defined as

Y c V c = ∂Y c ∂x Y c Ṽc = ∂Y c ∂x i α i Y i Ṽi , (18)
which also allows a Lewis number to be defined for the progress variable Le c = -λ/(ρ c p Y c Ṽc ). This information is essential for evaluating the diffusive flux of Y c , when using the premixed flamelet model in CFDsimulations. In this case, it was possible to store either Y c Ṽc or Le c as a function of Y c in a flamelet table. Using Eq. ( 17) instead of Eq. ( 8) in the flamelet transformation, the premixed flamelet equations for T , Y i and g c read in generalized form

g c c p ∂ ∂Y c g c λ ∂T ∂Y c + g c ∂ ∂Y c g c ρY c Ṽc ∂T ∂Y c + g 2 c k c p,k c p ρY k Ṽk ∂T ∂Y c -ωc ∂T ∂Y c + ωT c p = 0 , (19)
-g c ∂ ∂Y c g c ρY i Ṽi + g c ∂ ∂Y c g c ρY c Ṽc ∂Y i ∂Y c -ωc ∂Y i ∂Y c + ωi = 0 , (20)
-

g 2 c ∂ 2 ∂Y c 2 g c ρY c Ṽc -ωc ∂g c ∂Y c + g c ∂ ωc ∂Y c = 0 . ( 21
)
The consistency of these equations to the ones derived in previous works [START_REF] Nguyen | Combust[END_REF][6][7] is discussed in the supplementary material.

Boundary Conditions

The premixed flamelet equations are subject to the following boundary conditions

Y c = Y c,min : T = T 0 , Y i = Y i,0 , g c = 0 , Y c = Y c,eq : T = T eq , Y i = Y i,eq , g c = 0 ,
where the subscript 0 denotes the thermochemical state of the fresh mixture and eq marks the chemical equilibrium.

Numerical Methods

The premixed flamelet equations presented in the previous sections are solved with a hybrid Newton method [15] in a multi-grid approach with an in-house C++ flame solver [16]. After the equations are solved on a coarse grid, the grid is adaptively refined until a userspecified target. For the evaluation of diffusion coefficients and chemical source terms, the EGlib [START_REF] Ern | Multicomponent Transport Algorithms[END_REF][START_REF] Ern | EGLIB: a General-Purpose Fortran Library for Multicomponent Transport Property Evaluation[END_REF] and an in-house chemistry library are used, respectively. The GRI-Mech 3.0 [START_REF] Smith | [END_REF] is applied for methaneair flames and the mechanism by Varga et al. [20] is applied for hydrogen-air flames.

Results and Discussion

In this section, general characteristics of the premixed flamelet model are discussed first. Thereafter, flamelet results are compared to the corresponding physical space flame configurations for three premixed flame categories: C1 methane-air flames assuming equal diffusivities, In comparison to the physical space, the reaction zone is enlarged in progress variable space, while pre-heat zone and postoxidation zone are confined to thin layers at the boundaries.

C2 methane-air flames using a mixture-averaged diffusion model, C3 hydrogen-air flames using a mixture-averaged diffusion model.

The flamelet solutions are computed by solving Eqs. ( 19)-( 21) as previously outlined and the physical space solutions are obtained by solving Eqs. ( 1)-( 3) with the procedure described in [START_REF] Kee | PREMIX: a fortran program for modeling steady laminar onedimensional premixed flames[END_REF]. For all flames the equivalence ratio is varied between 0.8 and 1.2. Finally, the budget of the flamelet equation for the progress variable gradient is analyzed.

Premixed Flamelet Characteristics

In comparison to the representation of premixed flamelets in physical space, the numerical model in progress variable space exhibits different characteristics, as shown in Fig. 1. While the reaction zone is very thin in physical space, it almost spans the whole computational domain in progress variable space. Consequently, gradients can be accurately resolved in the most important part of the flame with fewer points compared to the physical space representation. In contrast to the reaction zone, the flame's pre-heat and post-oxidation zones are confined to thin layers at the boundaries in progress variable space. This poses numerical difficulties for an accurate and robust solution method for the premixed flamelet equations. In the pre-heat zone differential diffusion effects might occur while there is minimal reaction progress. To properly resolve the differential diffusion effects, especially for light species such as H 2 , the grid must be refined towards the pre-heat zone and the reaction progress variable should contain some of the light species. Otherwise, their profiles might exhibit large gradients in the pre-heat zone, which leads to errors in the laminar flame speed of above 5% for the

Y c = Y CO + Y CO 2 + Y H 2 O + 10 Y H 2 -Y O 2 .
flames studied in this work. In the post-oxidation zone, the gradient of the progress variable becomes small and slow chemical reactions control this region. Asymptotically, as g c → 0, diffusive terms vanish in the flamelet equations and the equations for the homogeneous reactor in progress variable space are recovered [21]:

∂Y i ∂Y c = ωi ωc and ∂T ∂Y c = 1 c p ωT ωc . ( 22
)
The above equation set allows only one boundary condition to be specified. Thus, even if g c is non-zero, but very small, the information about the post-oxidation zone (in the limit an equilibrium boundary condition) shows very little effect on the reaction zone. Alternatively to prescribing the (exact) equilibrium conditions in the post-oxidation zone, it is also possible to define the boundary condition as an open boundary, solving the premixed flamelet equations at this point. With this, the thermochemical state in the post-oxidation zone is a result of the computation as it is similarly the case for freely-propagating flames in physical space which recover close-to-equilibrium conditions at the flow outlet. Albeit this theoretical insight, both methods are interchangeable and for all computations discussed below the boundary conditions as outlined in Sec. 2.4 are used.

Temperature, Species and Gradient Profiles

Figure 2 shows the solutions of the temperature, the gradient g c and four species for the methane-air flames with equal diffusivities (C1). Overall, the premixed flamelet model with equal diffusivities accurately reproduces the physical space solution for all displayed variations of the equivalence ratio. For the cases discussed here, the freely-propagating flame solutions were computed with a multi-grid approach on 1000 points, while the progress variable space model was computed on 200 points (both adaptively refined from approximately 20 starting points). The progress variable is defined by

Y c = Y CO + Y CO 2 + Y H 2 O + 10 Y H 2 -Y O 2
for consistency with the computations with mixture-averaged diffusion (C2), which further accounts for differential diffusion. However, for the cases with equal diffusivities, a simpler definition of the progress variable would also yield similar results.

The results for the methane-air flames with non-unity Lewis numbers are displayed in Fig. 3. In accordance with the previous results, the premixed flamelet model accurately reproduces all solution quantities. Although they are very similar to the results from Fig. 2, differences due to differential diffusion are notable in the profiles of g c , CO 2 and especially H 2 . As mentioned earlier, the H 2 shows large gradients in the pre-heat zone due to differential diffusion in the pre-heat zone. If H 2 was not contained in the definition of Y c , it becomes very challenging for the premixed flamelet model to capture this profile accurately (not shown here). From this observation it becomes evident, that the accuracy of the model poses more constraints on the definition of the progress 0.8 0.9 1.0 1.1

1.2

Fig. 3: Comparison of premixed flamelet solutions (dot) to freely-propagating flame solutions (line) using a mixture-averaged diffusion model for methane-air flames for different equivalence ratios (0.8-1.2, blue). The progress variable is defined as

Y c = Y CO + Y CO 2 + Y H 2 O + 10 Y H 2 -Y O 2 .
variable beyond strict monotonicity. In Fig. 4 hydrogen-air flames are shown for different equivalence ratios. Although H 2 O is a common choice for the progress variable for these flames, again H 2 had to be included in its definition for a sufficient accuracy of the model. As already observed for the methaneair flames, the model correctly reproduces the physical space flame results. This is also true for light radical species such as H, which underlines the fact that the premixed flamelet model correctly captures differential diffusion effects.

The comparison of the freely-propagating flames to the premixed flamelet solutions illustrates, that the closure for the conditional strain rate Eq. ( 13) is exact and can be used to accurately predict the progress variable gradient g c . Inspecting the shape of all gradient profiles (C1-C3) shows, that there does not exist a generic profile for premixed flames, hence g c must be a solution quantity.

Comparison of the Laminar Flame Speeds

In progress variable space, the laminar flame speed s l can be computed by transforming Eq. ( 6)

s l = 1 ρ 0 - ∂ ∂Y c g c ρY c Ṽc + ωc g c . (23)
In this work, s l is evaluated at the normalized progress variable c = 0.5. A comparison of s l between the result from the progress variable space model and the physical space solution for all cases (C1-C3) is displayed in Fig. 5. The premixed flamelet model reproduces s l with less than 2.5% difference as compared to the physical space computations.

Analysis of the g c Budget

Figure 6 shows the budget analysis for the g c flamelet equation Eq. (21). It is observed that the diffusive term and the source term are important throughout the whole domain, while the Y c -drift term shows major influences only in the region with considerable reaction progress where c > 0.6. This characteristic coincides with the maximum of the progress variable source term, which is shifted to the right of the computational domain, analogously to the heat release rate in Fig. 1. The budget shows that both terms which contain the chemical

Y c = Y H 2 O -Y H 2 .
source term of the progress variable play a dominant role. Obviously, these terms would not be present if the control variable was a passive scalar, such as the mixture fraction for non-premixed flames. While nonpremixed flamelets require a forcing term in the form of strain to sustain a steady flame (also evident from the equation for the scalar dissipation rate, see [START_REF] Hasse | [END_REF]), this is not the case for premixed flamelets. However, it would make sense to also introduce such a term in the premixed flamelet model to extend it to capture influences of imposed strain. A corresponding method in physical space is the strong stretch theory for FGM-models [2].

Summary

In this work, a self-contained flamelet model for steady unstrained premixed flames is presented. Flamelet equations for temperature, species and the gradient of the progress variable are formulated and solved in progress variable space, which also yields a closure for the laminar flame speed. A generic transformation of the species diffusive flux is introduced, such that it can be computed with standard transport libraries on a Y c -grid. With this, the model can capture differential diffusion effects, since common diffusion models may also be used in progress variable space. The premixed flamelet model is then applied to three different flame categories, varying the equivalence ratio: methane-air flames assuming equal diffusivities (C1), methane-air flames using a mixture-averaged diffusion model (C2) and hydrogen-air flames using a mixture-averaged diffusion model (C3). It is shown that for all cases the model accurately captures the profiles of the solution quantities in accordance with the corresponding freelypropagating flame calculations in physical space. Furthermore, the progress variable space solution method can predict the laminar flame speed with differences less than 2.5% compared the physical space solution, where it is an eigenvalue.

The model presented here complements flamelet modelling approaches for the most important canonical premixed flame configuration. It remains a subject for future work to show how strain and curvature effects can be consistently incorporated into a progress variable space flamelet model and to explore how the model may be integrated into flamelet approaches for turbulent flames.

Fig. 1 :

 1 Fig.1: General characteristics of premixed flamelets -left: Temperature, progress variable and heat release rate (HRR) in physical space. Right: Temperature and heat release rate (HRR) in progress variable space. In comparison to the physical space, the reaction zone is enlarged in progress variable space, while pre-heat zone and postoxidation zone are confined to thin layers at the boundaries.

Fig. 2 :

 2 Fig. 2: Comparison of premixed flamelet solutions (dot) to freely-propagating flame solutions (line) assuming equal diffusivities (Le = 1) for methane-air flames for different equivalence ratios (0.8-1.2, blue). The progress variable is defined as Yc = Y CO + Y CO 2 + Y H 2 O + 10 Y H 2 -Y O 2 .

Fig. 5 :

 5 Fig. 5: Laminar flame speeds (C1-C3) for the physical space solution (lines with symbols) and for the premixed flamelet model (red dots).

Fig. 4 :

 4 Fig. 4: Comparison of premixed flamelet solutions (dot) to freely-propagating flame solutions (line) using a mixture-averaged diffusion model for hydrogen-air flames for different equivalence ratios (0.8-1.2, blue). The progress variable is defined as Y c = Y H 2 O -Y H 2 .

Fig. 6 :

 6 Fig. 6: Normalized budget of Eq. (21) over the normalized progress variable c for a stoichiometric methane-air flame. The balance represents the sum of the three terms.

Acknowledgments

Financial support is kindly acknowledged from the German Research Foundation (DFG) in the Project HA 4367/4-1.

List of Figures