Conformational Solvation Studies of LIGNOLs with Molecular Dynamics and Conductor-Like Screening Model

Molecular dynamics (MD) simulations were performed on sterically hindered α-conidendrin-based chiral 1,4-diols (LIGNOLs) from the naturally occurring lignan hydroxymatairesinol (HMR) using the GROMACS software. The aim of this study was to explore the conformational behaviour of the LIGNOLs in aqueous solution adopting the TIP4P model. The topologies of the LIGNOLs were constructed manually and they were modeled with the OPLS-AA force field implemented in GROMACS. The four most relevant torsional angles in the LIGNOLs were properly analyzed during the simulations. The determining property for the conformation preferred in aqueous solution was found to be the lowest energy in gas phase. The solvation effects on the LIGNOLs were also studied by quantum chemical calculations applying the COnductor-like Screening MOdel (COSMO). The hydration studies of the MD simulations showed that several of these LIGNOLs, produced from a renewable source, have a great potential of acting as chiral catalysts.


Introduction
The forest industries are developing new innovative products in addition to the traditional bulk products. A promising raw material for added-value products consists of lignans that are extracted in chemical pulping from residual knots.
The anticarcinogenic and antioxidative lignan hydroxymatairesinol (HMR) is found in large amounts in the knots of Norway spruce (Picea abies) [1]. It has been used in the synthesis of TADDOL-like α-conidendrin-based chiral 1,4-diols (LIGNOLs) [2] with the same functionality as TADDOLs [3,4] or BINOLs [5], which are often used as ligands for transition metal catalysed asymmetric synthesis. They have hindered structures containing two adjacent stereocenters, resulting in a fixed angle between the metal-complexing hydroxyl groups.
HMR has recently been studied by quantum chemical calculations [6,7] and by molecular dynamics (MD) [8]. The structures of the LIGNOLs included in this study have been quantum chemically optimized [2], and in this study MD simulations are used to explore the conformational changes of the LIGNOLs in aqueous solution. Such an analysis has to the best of our knowledge not been performed before.
Solvation effects are important to include in quantum chemical calculations because solvents play a role in the structures of many molecules. This can be done by using implicit solvation models such as, e.g., the COnductorlike Screening MOdel (COSMO) [9,10]. In a previous study [8] on the structurally quite similar lignan HMR, solvent effects calculated by COSMO were compared to solvent effects by the Polarized Continuum Model (PCM) [11][12][13], and COSMO was found to give more credible relative energies than PCM.
The three quantum chemically most stable conformers [2] of each of the LIGNOLs (three per stereoisomer for triphenyl) were chosen for this study. The four most relevant torsional angles in the LIGNOLs were properly analyzed during the simulations. Those are explained in detail in Figure 2, which also shows the numbering of atoms.

Solvation Effects in Ab Initio Calculations
It is well known that the surrounding medium might determine which conformation will be preferred. In order to corroborate this, the systems were studied in aqueous solution with the COSMO model at different levels of theory.
In Tables 1 and 2 the total electronic energies and the dipole moments, respectively, are expressed. The energies are given relative to the most stable conformer for each method used. The first two columns [2] show the total electronic energies or dipole moments in the gas phase. In the second column the level of theory used was B3LYP/TZVP. In the last two columns the values from the studies in the presence of water (ε r = 78.39) using the COSMO model are reported. The methods for the COSMO calculations in the third column is reoptimization at the restricted Hartree-Fock (RHF)/SV(P) level, and in the fourth column reoptimization with density functional theory (DFT) at the B3LYP/TZVP level of theory, i.e., the starting structures for the optimizations are the HF gas phase structures. The exact values of the torsional angles can be found in Table 3.   Table 3. The initial (HF) and the DFT optimized torsional angles.
Conformation Method α β γ δ As can be seen in Table 1 the energies from the COSMO calculations follow the same order as the conformers in the gas phase for diphenyl and tetraphenyl 1,4-diol. For triphenyl and tetramethyl 1,4-diol the case is different.
The trend, however, is the same that could be observed in the gas phase study in a previous work [2]. For the energetically more favourable conformers, a π − π interaction was formed between the phenyl ring at C7 and one of the phenyl rings at C9 . As stated in [2] the change was initialized by a turn of the phenyl ring at C9 for 3PhS3 and 3PhS7. In 3PhR3 the torsional angle α was turned backwards by ≈ 100 • , but the phenyl ring at C9 was to begin with in a correct position to form the desired π − π interaction. What is notable is that this trend can be seen already at the HF level (column 3) in the COSMO calculations.
For the tetramethyl conformer 4Met6 this could not be noticed at the HF/SV(P) level even if the DFT optimization changed the structure a little, though, causing a big energy yield. The phenyl ring at C7 tilts almost 50 • , then allowing a change in the aliphatic six-membered ring from an envelope like conformation to a boat conformation. This also forces one of the methyl groups closer to the phenyl ring at C7, exactly as it happened in the gas phase study [2].
When considering the dipole moments in Table 2, calculated with GAMESS at HF/6-31G* level, one can observe that the order between the conformers is preserved almost entirely for each method used. Generally the dipole moment always increases when a polar molecule is solvated. This is of course also true for the solvation models, as can be seen in Table 2.
When examining the values in Table 2 more carefully, one can see that the diphenyl 1,4-diols are the ones that follow the energy trend the closest. The only exception is µ DF T in gas phase for 2Ph9. For triphenyl and tetraphenyl 1,4-diol a trend can be seen that conformers with higher dipole moment are more stable than the other, especially according to the DFT calculations. This also holds for tetramethyl 1,4-diol, for which conformer 4Met6 has by far the highest dipole moment and is most stable in DFT.

Molecular Dynamics
The initial values of each of the four torsional angles mentioned in the Figure 2 are shown in Table 3. In Table 3 also the torsional angles in the DFT optimized structures [2] are shown, but the HF structures are the ones used as starting structures in the MD simulations. The torsional angles that change a lot during the DFT optimizations are marked in bold.      The first two simulations for diphenyl 1,4-diol were almost status quo. The four torsional angles stayed close to their starting values; α = 120 • or 300 • , β = 60 • , γ = 70 • and δ = 315 • , which can be seen in Table 3. In the first simulation a short conformational change occurs just after 3 ns into α = 150 • , δ = 225 • and γ < 60 • , and this occurs a few times between 7 and 9 ns.
The third simulation of diphenyl 1,4-diol, however, is perhaps the most interesting case of all in this study. The torsional angle γ is as usual the most stable one in these LIGNOLs, staying at the starting value 195 • . The highest populated conformation during this simulation is the one staying stable between 5 and 7 ns, with α = 120 • , β = 200 • and δ = 255 • . This value of δ was found to make it possible for the aliphatic six-membered ring to be in boat conformation, by that raising the stability of the tri-and tetraphenylated 1,4-diols in [2]. The starting value of β = 200 • , but this changes immediately to 300 • . The angle β changes back at 0.5 ns not causing any immediate change, but perhaps initializing the change of α → 150 • and δ → 255 • . At 1.5 ns α → 120 • and β back →< 200 • . At 2.5 ns δ → 300 • without any other consequences, but at 3 ns the conformation changes to α = 150 • , β = 300 • and δ ≥ 255 • . This does not stay for long until the conformation returns to α → 120 • and δ → 300 • , and soon afterwards, β back → 200 • . This pattern continues until the population of the conformation raises at 5 ns. After 7 ns δ changes again → 300 • , and at 8 ns the previous pattern from 3-5 ns starts again. If both α and δ change it seems to happen simultaneously, while β either initializes a conformational change or lags behind.

Triphenyl(R)
The triphenyl 1,4-diol seems to be a much more hindered molecule than the diphenyl, as it also should be, especially the R stereoisomer which is discussed here. In all three simulations δ stays most of the time at ≈ 255 • -the stabilizing value stated in [2]. The torsional angle α has the possibilities of 150 • or 330 • , with an exception in the second simulation, while α → 285 • and δ → 300 • at 9 ns, also raising β and γ with a few degrees, and the same quickly back and forth at 6.5 ns.

Triphenyl(S)
In the first two simulations of the triphenyl(S) stereoisomer, the conformations stay quite unchanged at α = 150/165 • , β = γ = 195 • and δ = 255 • . Just before 5 ns in the second simulation α → 105 • , but this does not seem to affect the rest of the torsional angles studied.

Tetramethyl
A bit surprisingly tetramethyl 1,4-diol is the most hindered molecule of the LIGNOLs that were studied. Even the conformations are quite similar in the three simulations: α = 150 • , β = 180/300 • , γ = 75/195 • and δ ≥ 240 • . In the first simulation, though, α changes three times → 105 • (at 1 and 3.5 ns), but this has very little effect on the rest of the studied angles.

Tetraphenyl
The originally planned end product in the synthesis route of the LIGNOLs, i.e., tetraphenyl 1,4-diol, is a more interesting case from the torsional point of view than the last one, even though, β = γ = 195 • through all of the simulations. However, α and δ fluctuate most of the time between the conformers α = 135/330 • −δ = 255 • , and α = 105/285 • −δ = 285/300 • . In the second simulation just before 7 ns, α flips drastically → 135 • , though not causing anything at all on the rest of the studied torsional angles.

Discussion on the Simulations
For diphenyl 1,4-diol, the highest populated conformations seemed to be 2Ph1 and 2Ph2, as can be seen in Table 2 in [8] and in Table 3 in this article. Another quite highly populated conformation in the simulations in TIP4P water is 2Ph9, which is the third one picked out as starting geometry for the simulations, but with the exception that the torsional angle δ is rather the favourable 255 • than the starting value ≈ 315 • . However, this change of δ seems to imply that α takes the uncommon value of 150 • , also causing or implying β to be 300 • .
In the case of triphenyl 1,4-diol, especially the R stereoisomer, it is more comprehensible that the value 150/330 • occurs for α as there is an additional phenyl ring substituted on C9. The highest populated conformers could be said to be ≈ the DFT optimized 3PhR3, and 3PhS3 and 3PhS7, respectively, even though the values for the S stereoisomers vary a bit from the optimized ones. The small variation in β and γ from the optimized ones might, though, be explained by the sensitivity of them, when forming hydrogen bonds between the hydroxyl groups, O9-H and O9 -H, and the TIP4P water. The torsional angles in the highest populated conformer in the last triphenyl(S) simulation is actually very close to 3PhS10.
In tetramethyl 1,4-diol it is simply the three most stable conformers in gas phase that dominate, i.e., 4Met2, 4Met3 and 4Met6 , the last two of which are actually the same conformer.
The last molecule studied, tetraphenyl 1,4-diol, is trickier to summarize in the system of the conformers optimized in gas phase. The most frequent value of δ in the simulations, 255 • , only occurs for the conformers 4Ph5 and 4Ph6, but for those the other torsional angles vary a lot. The only conformer of the simulations, fitting well with the gas phase values, is 4Ph4 and 4Ph8, actually the same, which occurs in the second simulation. The big fluctuations compared to the other LIGNOLs might, though, be a sign of various populations of the molecule.
Opposite to the conclusion for HMR [8], the quantum chemical calculations (in gas phase), thus, seemed to predict reliable trends for how these molecules act in solvents. Perhaps this was more reliable as the molecules were bigger than HMR.

Hydration Effects
In the LIGNOLs there are four hydrogen bonding acceptor oxygen atoms in methoxy groups. However, the most interesting oxygens from the reaction point of view are those in the metal-complexing hydroxyl groups, O9-H and O9 -H. In those groups there are also two hydrogen bonding donors.
In order to understand the hydration effect more properly, the g_hbond analysing program implemented in GROMACS was used to study the number of hydrogen bonds for the oxygen atoms O9 and O9 , and totally for each LIGNOL conformer, as well as the average lifetime of the uninterrupted hydrogen bonds, which are all shown in Table 4. In a previous study [8] the average number of hydrogen bonds per time frame for the oxygen atoms in the methoxy groups was ≈ 0.24, which means that the four methoxy oxygens in the LIGNOLs contribute with one hydrogen bond altogether. The rest consists mainly of contributions from the reactive centre, i.e., the hydroxyl groups, O9-H and O9 -H.
The first three columns in Table 4 show that O9 seems to be a bit more likely to form hydrogen bonds to the solvent. This conclusion might, though, be erroneous due to the fact that the hydroxyl groups O9-H and O9 -H form an internal hydrogen bond, which takes away one connection to the solvent, for those conformers that have the OH groups pointing in the same direction. However, it is notable that tetramethyl 1,4-diol is more likely to form hydrogen bonds to TIP4P and tetraphenyl is less likely so, mainly due to the small tendency of O9 to form hydrogen bonds to TIP4P.
Considering the lifetimes, in [8] it was stated that the average lifetime of the uninterrupted hydrogen bonds was calculated to be 1.2 ps. The last column shows that this also holds quite well for the LIGNOLs. A few are just below 1 ps. When looking at columns 4 and 5, more interesting values can be noticed, especially for O9 (col. 5). A correlation can, however, be seen to the number of hydrogen bonds as the lifetimes are longer for tetramethyl 1,4-diol and shorter for tetraphenyl. A shorter lifetime for a large average number of hydrogen bonds may imply that they are quite weak, meaning that the hydrogen bonds from O9 in 2Ph1 and 2Ph2 might be strong. This again could be very important for the application of these LIGNOLs as metal-binding agents, as the bonding to a metal-atom catalyst would act as the hydrogen bonding to TIP4P water. Diphenyl 1,4-diol is the only LIGNOL in this study with phenyls at C9 and not at C9, so the reason for this phenomenon is probably the electronic effects of the phenyl rings at C9 .
The hydrogen bond lengths were also analyzed by using g_hbond, and the mean value of the hydrogen bond lengths is approximately 0.28 nm, exactly as in [8].
The MD simulations were performed using GROMACS version 4.5.3 software [30][31][32][33][34]. Water was described using the TIP4P model [35], and the LIGNOLs were modeled with the OPLS-AA force field [36] implemented in GROMACS. The topologies of the LIGNOLs were constructed manually, and they comprised 415 (2Ph), 474 (3Ph), 533 (4Ph) and 369 (4Met) internal coordinates, respectively. In order to get reasonable atomic charges to help for choosing suitable atom types with the hand-tuned charges available in the force field, electrostatic potential fit (ESP) charges were studied with GAMESS at HF/6-31G* level. For O9 and O9 (shown in Figure 2) the OPLS atom type opls_154 with the atomic charge −0.683 was found to be the most suitable one, and for the other four oxygens (O3, O4, O4 and O5 ) the atom type opls_179 with the atomic charge −0.285 was chosen. An important detail to consider is also that the sum of the atomic charges in a charge group should be an integer or equal to zero.
The initial (quantum chemically optimized) conformations of the LIGNOLs were taken from our previous work [2]. Each conformation was placed at the center of a cubic box with the dimension between 5.2 and 5.6 nm (volume = 144-174 nm 3 ) and solvated by 4802-5795 water molecules. Each system was first energy minimized < 2000 kJ mol −1 nm −1 using steepest descent for 3-121 steps. Then the system was equilibrated at 398 K for 50 ps, and finally the production simulation was run for 10 ns with the temperature maintained at 298 K using the Berendsen thermostat [37]. The pressure was maintained at 1 atm using the Berendsen barostat [37]. A 1 fs time step was used in all simulations.
A cutoff of 0.9 nm was applied to short-range nonbonded interactions, and for long-range electrostatic interactions the particle mesh Ewald (PME) method [38,39] was used with grid spacing of 0.12 nm and fourth-order interpolation. In all simulations system snapshots were collected every 500 steps, i.e., 0.5 ps, for subsequent analysis. In this time only electronic excitations and bonding vibrations will occur, but those can be ignored when studying the conformational preferences of the system.

Conclusions
• In MD simulations on the LIGNOLs, the conformations preferred were the energetically most favourable ones according to quantum chemical DFT calculations in gas phase, almost irrespective of the dipole moment. • The four most relevant torsional angles α − δ, defined in Figure 2, varied during dynamics in accordance with their symmetry. • The torsional angle δ was generally more preferred at the stabilizing value 255 • than what was seen in the gas phase optimizations in [2]. • No strong correlation patterns were found, but in the last simulation of 2Ph9, α and δ changed simultaneously, while β either initialized a conformational change or lagged behind. • In the hydration studies 2Ph1 and 2Ph2 were found to have strong hydrogen bonds from O9 , which could be very important for the application of these LIGNOLs as metal-binding agents. • The hydration studies of the MD simulations show that several of these LIGNOLs, produced from a renewable source, have a great potential of acting as chiral catalysts.