Frequency Dependent Non- Thermal Effects of Oscillating Electric Fields in the Microwave Region on the Properties of a Solvated Lysozyme System: A Molecular Dynamics Study

The use of microwaves in every day’s applications raises issues regarding the non thermal biological effects of microwaves. In this work we employ molecular dynamics simulations to advance further the dielectric studies of protein solutions in the case of lysozyme, taking into consideration possible frequency dependent changes in the structural and dynamic properties of the system upon application of electric field in the microwave region. The obtained dielectric spectra are identical with those derived in our previous work using the Fröhlich-Kirkwood approach in the framework of the linear response theory. Noticeable structural changes in the protein have been observed only at frequencies near its absorption maximum. Concerning Cα position fluctuations, different frequencies affected different regions of the protein sequence. Furthermore, the influence of the field on the kinetics of protein-water as well as on the water-water hydrogen bonds in the first hydration shell has been studied; an extension of the Luzar-Chandler kinetic model was deemed necessary for a better fit of the applied field results and for the estimation of more accurate hydrogen bond lifetime values.


Introduction
The dielectric properties of protein solutions have been the object of several studies over the years [1][2][3][4][5][6][7][8][9][10]. Although an increasing amount of knowledge has been accumulated regarding their dielectric response in a broad frequency range, the wide application of microwave technology in the telecommunication branch and in medical diagnosis and therapy [11][12][13] has made the study of the biological non thermal effects, specifically in the microwave band, of particular interest. The dielectric spectra of a protein can be understood generally in terms of PLOS

Theoretical Background Dielectric function
Here we assume that the applied oscillating electric field is parallel to the z-axis and is homogeneous in the entire space occupied by the simulation cell. Commonly, the oscillating field is given by: (Symbols have their usual meaning) Now let us consider a system of m different components, each of them (k = 1,. . ., m) behaving as a dielectric material. Their induced dipole moments at frequency ω can be shown to follow the expression [25]: M k ðo; tÞ ¼ A k ðoÞ Á cosðotÞ þ B k ðoÞ Á sinðotÞ ð2Þ The frequency dependent constants A k and B k are connected directly with the real and the imaginary part of the dielectric constant [25], via equations: Here V is the simulation cell volume. The total dielectric constant of the system is simply given by the summation of all the above m dielectric terms, that is: It is important to mention that in this work we do not include the dielectric behavior of ions. This is actually a common practice in dielectric relaxation spectroscopy experiments when a better resolution of the β, γ and δ-processes is required. This is usually realized by removing the DC conductivity of ions [6], i.e. by subtracting the term iσ 0 /E 0 ω from the dielectric spectra. Moreover, as we showed in a previous work [14], due to the small number of ions, their contribution to the dielectric constant is only minor and can be safely disregarded.

Hydrogen bond dynamics
The dynamics of hydrogen bond (HB) formation between water and the protein residues as well among water molecules themselves, can be well characterized by two time correlation functions [39], the continuous hydrogen bond time correlation function S(t) and the intermittent hydrogen bond time correlation function C(t). These correlations functions are defined by the following equations: Where h(t) and H(t) are two hydrogen bond population variables. Specifically h(t) is unity when a particular pair of sites is hydrogen bonded at time t and zero otherwise. On the other hand H(t) is unit when a specific pair remains continuously bonded from time 0 to time t. Thus S(t) describes the probability that a hydrogen bond has formed and remains stable at all times up to t (a strict definition of lifetime), while C(t) describes the probability that a tagged pair is hydrogen bonded at time t, given it was bonded at time zero, i.e. C(t) allows breaking and reformation of hydrogen bonds of a tagged pair at intermediate times.
Appropriate multi exponential fits on the above dynamic functions lead to the calculation of the so called intermittent τ C and continuous τ S average relaxation times of the proteinwater or inter-water system. These time averages provide a good description of rotational and translational relaxation motions as well as of kinetic bond mechanisms in different parts of the system (especially in active sites and hydration shells). Following the concepts and notation of reference [39] the kinetics of breaking and formation of protein-water bonds can be described as: Where B and QF, are the bound and quasi-free state of a specific pair of sites that participate to form a hydrogen bond. Since QF is defined by a cutoff distance, in our implementation water molecules arrive to this state either by rotating around the protein site or around their center of mass. A simple rate equation for the "reactive flux" has the following form [40,41,42]: Where k 1 and k 2 are the forward (breaking) and backward (re-formation) rate constants for the bound and quasi-free state respectively. C PW (t) is the intermittent hydrogen bond time correlation function C(t) in case of protein-water HBs. N pw (t) is a time correlation function that takes into consideration the diffusion of water molecules and can be calculated by the formula: Where H'(t) is unity when a pair of sites is closer than a cutoff distance R PW (usually taken equal to 3.2 Å) for protein-water interactions at time t and zero otherwise. I.e., the N pw (t) correlation accounts for water populations which are not any more hydrogen bonded to protein but remain close enough to reform a broken bond. The above model is first introduced by Luzar and Chandler (L&C) [40][41][42].
Extension of luzar and chandler model of hydrogen bonds kinetics. In this section we modify Eq (6) in order to include an additional relaxation process. This modification is necessary to describe better our simulation results related to the influence of the electric field on the hydrogen bonds kinetics. The calculation of the dynamic functions starts considering the instant configuration of the hydrogen bond network between protein and its neighboring water molecules, which defines the t = 0. For the above system one can write the following double reversible reaction: Following the fate of each previously bound water molecule, the left part of reaction (9) describes the direct transition from B to a second quasi free state QF' and vice versa (Fig 1). That is, we consider the possibility of some originally bound water molecules to be found slightly beyond the cut-off distance R PW from a specific protein site after the initial breaking. On the other hand QF' represents a second "activated state" along a reaction coordinate defined by the distance (very close to R PW ) from the protein site allowing some water molecules to escape, while others to directly reform the HB. In our calculations we consider a limit for this distance determined by the largest distance R max observed after the first bond breaking, which is actually related to the frame recording rate. All previously bound water molecules found in a spherical shell between R PW and R max belong to the small QF' population not captured by the L&C model. The right side of reaction (9) refers to the transition between the bound B and quasi-free state and vice versa. According to model (9) a "reactive flux equation" can be written for C PW similar to Eq (7): where k 1 accounts for the total (k 1 +k' 1 ) forward (breaking) of HBs of model (9). F PW (t) is a function introduced to describe the time correlation of QF' water molecules. The F(t) correlation is given by: Where h(t) and H ' (t) are the previously defined hydrogen bond population variables. Function Q takes at time t the value 1 if there is at least one water molecule (after an initial HB breaking) at distance R PW < r R max (t 0 ) from a protein's site able to reform directly the previous broken hydrogen bond. Here t 0 defines the starting simulation frame for the correlation analysis.
Since we count hydrogen bonds using labeled water molecules from concrete populations, their numbers will become gradually less and less resulting to a relaxation of the corresponding HBs correlation functions C(t), N(t) and F(t). According to our findings, the first hydration shell is less dense under field conditions than without. Therefore, water molecules in the QF' state have a smaller probability to defuse away than under zero field conditions allowing for a higher HBs reformation probability. It is this behavior which we incorporated in the reactive flux equation. As it will be shown in section "Results and Discussion" the incorporation of the F dynamics into the basic kinetic formula of L&C leads to better least square fits and consequently to more accurate calculations of average lifetimes especially in the case of electric field induced effects. In this work k(t) is calculated either as -ΔC/Δt from the simulation data giving k sim , or from a regression fit of two independent variables (C, N) adjusting the rate constants k 1 , k 2 to satisfy Eq (7), resulting to k LC . When determined from a three variables (C, N, F) fit to satisfy Eq (10) k(t) is named k F . From the rate constant k 1 determined by the fits to Eqs (7) or (10) we can then calculate relaxation times z = 1/k 1 for protein-water (PW) and inter-water (WS1) hydrogen bond breaking.

Simulation Details
A series of MDS has been performed aiming to a) the determination of the dielectric spectrum of a lysozyme solution with a direct method (application of an oscillating electric field at different frequencies of the spectrum) and b) to the study of possible non thermal effects of microwave radiation on lysozyme's structural and dynamical features.
As reported in our previous work [14], the starting complete structure of lysozyme (PDB entry: 2LZT) has been solvated in an orthogonal box (48.10×52.40×66.65 Å 3 ) with 4672 TIP3P water molecules giving a protein concentration of 9.88 mM (w/w = 14.5%, mass of H 2 O over dry protein = 5.88). Disulphide bonds have been added at Cys6-Cys27, Cys30-Cys115, Cys64-Cys80 and Cys76-Cys94. The positive charge of lysozyme (+8) has been neutralized by the addition of ions (18 Na + and 26 Cl -). The molecular dynamics program NAMD [43] and the CHARMM force field [44] have been used to energy minimize and pre-equilibrate the system for 2.5 ns followed by a production run of 30 ns with Langevin thermostat at 300 K in the NVT ensemble. The SHAKE algorithm [45] has been applied to all bonds and interactions were calculated with a cut-off radius of 16 Å. Electrostatic interactions were calculated with the Particle Mesh Ewald technique [46] while assuming tin-foil boundary conditions. The time step of the simulations was 2 fs and configurations of the system were stored every 0.5 ps.
Starting from the fully equilibrated state the above system was exposed during simulation to a uniform internal cosinusoidal electric field (E 0 = 4.2 mV/Å) at 21 frequencies in the range of 10 7 to 10 11 Hz using the eField parameters in the NAMD configuration file [47]. For the sake of comparison, the applied field is in the same order of magnitude as the field across a biological membrane (~1 mV/Å), while it is five orders of magnitude higher than that emitted by a cell phone (~10 -5 mV/Å). However such high intensities are necessary for small effects to exceed the thermal noise within the limited time scales generally amenable to molecular simulations. Field effects, like those presented later in this study, may arise as result of much lower intensities after a long time (minutes or hours) exposition. Unfortunately, such time domains are currently inaccessible to MDS. The simulations were performed in the NVT ensemble at 300 K. The computational time of each simulation was frequency dependent, i.e. for the lowest frequency of 10 7 Hz a total simulation time of almost 100 ns was necessary in order to cover a full cycle of the field (data saved every 2 ps). With increasing frequency the simulation time was decreased reaching a value of 0.01 ns at 10 11 Hz (data saved every 20 fs). For each frequency above 10 9 Hz twenty independent simulations were carried out in order to improve the statistics for our calculations.
Prior to applying a frequency to the protein solution the highest field strength value has been determined for which the system responds linearly. This is necessary in order to avoid break down of the linear response theory. A series of simulations have been performed with a static electric field applied on a system consisting of water molecules. The field amplitude has been increased starting from a value of 0.42 mV/Å and the upper acceptable limit for which the total dipole moment of water responds linearly has been determined to be 4.2 mV/Å. This value has been adopted for further frequency dependent simulations.
The simulation data have been analyzed using home made Tcl scripts running on the VMD program [48]. Further, the Gnuplot program has been used for fitting of the respective curves.

Results and Discussion
Dielectric calculations based on application of oscillating electric fields The dipole moment of the system components has been calculated for each one frame of the simulation trajectories derived for different applied frequencies of the dielectric spectrum. Though a field amplitude of 4.2 mV/Å is low enough to guaranty linearity, the price to pay is a low signal to noise ratio in the case of a protein. It is apparent from Fig 2A that the z-component of the dipole moment response of the total water (red line) at 12.5 GHz is periodic with time, giving a very good signal compared to noise. A fit to Eq 2 (blue line), provided the values for A k (ω) and B k (ω) which were then introduced in Eqs (3a) and (3b) to calculate the real and imaginary part of the dielectric function at this frequency. The fact that 12.5 GHz is at the absorption maximum of water, but far enough from that of protein, leads to a different behavior of the protein ( Fig 2B). Apparently, in the case of only one simulation at this frequency (green line) the statistic is too poor to reveal a clear signal, since thermal motion of the surface residues overcomes the response to an electric field of the given strength. For this reason 20 independent simulations have been averaged (red line) reducing drastically the noise. This allowed an acceptable fit (blue line) to the analytic function of Eq (2). It should be noted that the number of simulations needed to eliminate noise depended on the frequency and decreased gradually with decreasing frequency. Only one simulation is performed at 10 7 Hz, with a total simulation time of 100 ns for one full field cycle. This procedure spanned the dielectric function of the system and its components over a frequency range of 10 7 -10 11 Hz (Figs 3 and 4).
The calculated dielectric functions of a lysozyme solution and of protein alone are illustrated in Figs 3 and 4 respectively. Blue lines represent the dielectric function as calculated by the F-K approach and LRT according to a detailed decomposition method [14], while red squares represent values of the frequency dependent real and imaginary parts of the dielectric constant, as calculated by applying an oscillating electric field. In the case of protein alone the chromatic code is green lines and orange triangles (see Fig 4). Both, the locations and the heights of the peaks practically coincide. Moreover, concerning the δ-relaxation process (10 8 −10 9 Hz) the agreement is very good too, supporting the bi-modal nature of its dispersion  [14]. The real (a) and imaginary (b) part of the dielectric function of the entire system as calculated by the F-K approach (blue line) as well as by applying an oscillating field (red squares).
doi:10.1371/journal.pone.0169505.g003 curve [6,14]. Our calculations of the dielectric function by applying oscillating fields shows that the spherical protein approximation in the Fröhlich-Kirkwood approach works very well in the case of an ellipsoidal protein like lysozyme and keeping the field intensity lower than the saturation limit. Differences in the calculated dielectric spectra were expected in studies of proteins of shapes deviating rather largely from the spherical symmetry. However, in order to confirm the later statement, more studies with different proteins have to be performed. This is a key result, because it demonstrates the ability to calculate the dielectric properties of proteins from atomistic simulations without the assumption of an almost spherical symmetry assumed by the F-K approach. Moreover, the versatility of this non-equilibrium approach (at least for field strengths as high as the one used here) allows the study of the field effects on various protein structural and dynamic properties, as will be discussed in the next section.

Field effects on protein structure
In this section the effect of the applied field on protein's structure will be examined in detail. Protein structural changes will be considered, as reflected in the time course of protein's RMSD (C α root mean square deviation relative to lysozyme's crystal structure), in the RMSF  [14]. The real (a) and imaginary (b) part of the dielectric function of the protein as calculated by the F-K approach (green line) as well as by applying an oscillating field (orange triangles). (root mean square fluctuations) of individual residues and in the % of protein's secondary structure. Fig 5 displays the time evolution of protein's backbone deviation from the crystal structure at four frequencies (f 1 = 10 7 Hz, f 2 = 3×10 8 Hz, f 3 = 3×10 9 Hz and f 4 = 2×10 10 Hz) and at zero field. Frequency f 1 of the main protein absorption as well as f 2 next to it, exhibit protein structure deviations from the crystal structure relatively early (after 20 ns), while f 3 after about 40 ns. On the other hand, frequency f 4 of the bulk water absorption, well beyond the main absorption maximum of protein, shows a delayed (after 70 ns) but significant structural deviation larger than that appeared at f 3 and similar to that characterizing f 1 . According to references [37] and [38] the application of more intense oscillating fields (50 or 500 mV/Å) can lead to a dramatic increase in the RMSD values at frequencies away from the β-process zone.
In relation to secondary structure (Fig 6), VMD recognizes in the crystal structure six helical segments, helix-1: Arg5-His15, helix-2: Tyr20-Gly22, helix-3: Leu25-Ser36, helix-4: Cys80-Leu83, helix-5: Ile88-Asp101 and helix-6: Val109-Cys115 as well as five beta segments, strand-1: Thr43-Arg45, strand-2: Thr51-Tyr53, strand-3: Ile58-Asn59, strand-4: Cys64-Asn65 and strand-5: Ile88-Pro89. At the end of the simulation without field the same software recognizes helix-1, helix-3, helix-4, helix-5 and helix-6 and three beta conformations, strand-1, strand-2 and strand-3. A comparison of protein's secondary structure before and after exposure to the field at different frequencies (Fig 6) reveals changes in the amino acid ranges both of helices and β-strands without affecting significantly the overall structural organization of lysozyme. Specifically, all simulations lack short helix-2 and β-strands -4 and -5. On the other hand, it seems that field conditions promote the formation of helix-7 which is absent in the crystal structure and without field. Moreover, strand-3 is missing at 10 7 Hz, while strand-2 is missing at 3×10 9 and 2×10 10 Hz. Some elongation and shortening effects can also be observed. Frequency f 1 destabilizes (Fig 7) the right end of helix-3 including Glu35 and the left end of strand-2 including Asp52, both being the key catalytic residues of lysozyme, as well both ends of helix-5 (Ile88-Asp101). Table 1 summarizes the alterations in secondary structure content-beside radius of gyration (R g ) and RMSD. Since the numbers given in this table are averages over 5000 values the standard errors of the mean are very small (< 0.034). It is obvious that the field promoted the helical conformation in the C-terminus. The above changes in protein's secondary structure had no effect on the radii of gyration and do not go hand in hand with effects on RMSD. Besides, the largest RMSD is observed at 2×10 10 Hz, which is the main absorption frequency of bulk water. More substantial structural changes may be expected in proteins with less disulphide bonds.
Considering the above information, it is interesting to see which protein regions experience the most significant structural fluctuations under the influence of the oscillating field and if  Table 1. Alterations in percentage of secondary structure content, radius of gyration (R g ) and RMSD of lysozyme due to application of an oscillating field at 10 7 Hz, 3×10 8 Hz, 3×10 9 Hz and 2×10 10 Hz compared to zero field conditions and the crystal structure. All parameters are calculated as averages over the last 1ns (5000 frames, every 0.2 ps) of the corresponding 100 ns simulations. *RMSDs are calculated relative to the crystal structure of lysozyme (PDB:2LZT).

Non-Thermal Effects of Microwaves on a Lysozyme Solution: A Molecular Dynamics Study
there is a frequency effect. This is shown in Fig 7, which depicts the excess C α -atom fluctuations ΔRMSF (after subtracting the zero field RMSFs) at frequencies f 1 , f 2 , f 3 and f 4 . Interestingly, different frequencies affect different protein segments and the most striking fluctuations are located in loop regions, except one on the right end of helix-2 (Lys25-Ser36) triggered by f 1 and which contains Glu35, an important residue for the catalytic function of lysozyme. This observation differentiates from findings corresponding to more intense fields [38], where RMSF curves of 10 mV/Å at 2.45 GHz are found to be almost identical relative to zero field conditions. On the other hand, in case of an even higher and static field (50mV/Å) two charged segments (Thr43-Ile55 and Asn65-Leu75) showed high perturbations and disruption of hydrogen bonding, while segment Leu75 to Leu129 was not affected. According to this picture, the binding affinity of lysozyme to its substrate could be affected by the field at frequency f 1 . An indication of such an effect can be distinguished in the enhanced mobility of Phe34 (ΔRMSF % 2 Å) and Leu56 (ΔRMSF % 2.7 Å) which belong to the substrate's binding cleft. The structural effects described above combined with effects of the field on the intra-protein and protein-water hydrogen bonds -next in the text-may be correlated to effects of long exposure to cell phone radiation on enzyme activity, as observed in the case of acetylcholinesterase [49].

Hydrogen bond analysis
Next, the number and the dynamics of intra-protein (PP), protein-water (PW) and inter-water (WS1) hydrogen bonds in the first hydration shell (S1) will be discussed in detail. Here the realization of HBs has to satisfy the following geometric criteria: a cutoff distance of 3.2 Å for PW HBs and 3.5 Å for WS1 and PP HBs, while the angle formed by the donor, hydrogen, and acceptor is within 150-180 degrees.
Hydrogen bond dynamics. The number of hydrogen bonds in different subsystems of the lysozyme solution does not seem to be affected dramatically by the oscillating electric field at any frequency studied here (Fig 8). However, the small differences observed (Table 2) between the average number of hydrogen bonds over the last 10 ns of the 100 ns trajectories can be considered statistically significant, because of the relatively large number of the samples used (5000 values each). The most notable effect is observed in frequency f 1 (10 7 Hz), where it appears that the number of intra-protein hydrogen bonds decreased (14.3%) relative to zero field conditions. An opposite effect (increase) is shown for higher frequencies, with the largest (13.1%) observed at f 2 (3×10 8 Hz). Eventually, field effects on hydrogen bonding patterns become more prominent if one considers the number fluctuations of HBs (Fig 9) in conjunction with their dynamics (Fig 10, Tables 2-4). Fluctuations are calculated as time dependent variances of the number of HBs. A cursory glance of Fig 9A reveals a strong effect of frequencies f 1 and f 3 on the intra-protein HBs fluctuations, while fluctuations at f 2 and f 4 do not deviate from zero field conditions. Hydration water plays a significant role in all protein interactions with other molecules of biological interest. Therefore, it is very interesting to examine the behavior of HBs between protein and its surrounding water molecules (PW) in the presence of the field. A visual inspection of Fig 8B shows systematic differences between the number of protein-water HBs at different field conditions calculated for the last 10 ns of the 100 ns trajectories. Accordingly, Table 2 shows statistically significant decreases relative to zero field at f 2 , f 3 (δ-relaxation) and f 4 (γ-relaxation), with the largest (4%) observed at f 4 . For the same reason as in Table 1 the standard errors of the mean are once again very small (<0.24). Moreover, inspecting the protein-water HBs fluctuations (Fig 9B), it seems that at frequencies f 2 and f 4 the deviations from the mean converge at earlier times and to lower values compared to the behavior describing f 1 and f 3 and in the absence of the field.
Differences are more clear in WS1 inter-water hydrogen bonds. There is a small but distinguishable trend of the field to decrease the number of inter-water HBs (Table 2), with the largest decrease observed at f 4 (9.6%). Since PW-HBs are formed between protein and water molecules belonging to S1, it might be expected that inter-water HBs fluctuations in S1 would be related to those of PW-HBs (Fig 9B and 9C). However, they should not be same because S1 contains also water molecules which may not be hydrogen bonded to protein. Noticeably, the number of HBs (PW+WS1) per water molecule in S1 remains equal to 1.07 at all field conditions.
Another result related to the above picture concerns the impact of the field on the density of the first hydration layer. The volume of S1 water layer required for these calculations have been determined by approximating lysozyme as an ellipsoid [50]. A thickness of 3 Å has been assumed for S1, as determined by combined SAXS and SANS experiments [51]. Our findings suggest a higher density (13.6%) of the first hydration layer than that of bulk water (assumed to possess a density of 996.516 kg/m 3 at 27˚C), very close to experimental results (10%) derived using SAS [51]. Interestingly, we found that the presence of the field imparts a decrease of the excess density, with the lowest one observed at 2×10 10 Hz (Table 5). This correlates with the lower number of PW and WS1 HBs at this frequency ( Table 2).
The above frequency dependent changes may be explained more likely in terms of an interplay between two hypothetical mechanisms, direct field induced HB breaking and HB formation affected by field facilitated diffusion. Given that the density of S1 is normally higher than those in S2 and bulk, a higher water mobility induced by the field (lower friction) can result to a diffusion of S1 water molecules into the S2 or the bulk region. This would reduce the number and density of S1 water molecules, in accordance to our observations (Fig 8D, Table 5). A smaller number of S1 water molecules would reduce the probability of HBs formation between themselves and between protein and water ( Table 2). As a consequence (at high frequencies) some of the charged and polar amino acids, which have lost their water partners, will form ion bridges and therefore PP HBs are favored (Table 2). On the other hand, at lower frequencies, the field seems to affect PP HBs more directly by reducing their number and thereby increasing PW HBs (Fig 8A and 8B, Table 2). The hypothesized higher mobility induced by the field is probably connected to a momentary increase in temperature, which is corrected by the Fig 8. Field effects on the number of hydrogen bonds and S1 water molecules. Number of hydrogen bonds and S1 water molecules during the last 10 ns of 100 ns simulations in oscillating electric fields at f 1 = 10 7 Hz (β-process), f 2 = 3×10 8 Hz, f 3 = 3×10 9 Hz (δ 2 -process), f 4 = 2×10 10 Hz (γ-process) and without field. (a) Intra-protein HBs. Frequency f 1 has a clear negative effect on the number of intra-protein HBs, while f 2 and f 4 a small positive one. (b) Protein-water HBs. The oscillating electric fields reduce the number of protein-water HBs. (c) Water-water HBs in S1. The field decreases the number of inter-water HBs. (d) Number of S1 water molecules. Frequencies f 2 , f 3 and f 4 decrease the number of water molecules.
doi:10.1371/journal.pone.0169505.g008 Averages and standard deviations (SD) are calculated over the last 10 ns of the 100 ns trajectories (5000 values, every 2 ps). HB life times are calculated as ζ = 1/k 1 . The statistical significance has been calculated for differences between the average values of Table 2 and those under zero field conditions. All differences are found to be statistically significant (P-Values < 0.0001) at a significance level < α = 0.001.
doi:10.1371/journal.pone.0169505.t002 thermostat in the next steps. Of course this hypothesis has to be verified and expressed in mathematical terms, which will be the task of next studies."

Kinetics of hydrogen bond breaking and formation
The last part of this section deals with the dependence of the dynamical properties of system's hydrogen bonds on frequency. In recent years several simulation studies focused on the dynamics and kinetics of hydrogen bonds based on calculations of various correlations functions [52]. Concerning the first solvation shell S1 and various active sites of proteins, a different relaxation behavior of PW hydrogen bonds has been found due to the structural heterogeneity of the different protein segments [52]. In addition, in water simulations with application of a static electric field, it has been found that the kinetics of HBs are affected for field strengths above 100 mV/Å in agreement with theoretical findings, suggesting that the average number of hydrogen bonds remains constant up to the aforementioned field limit [53]. In the following paragraphs we analyze HBs kinetics at frequencies f 1 , f 2 , f 3 and f 4 and at zero field conditions in two hydrogen bond systems, protein-water and inter-water in S1. Protein-water hydrogen bonds. In order to make a comparison between the kinetics of bond breaking and reformation concerning the protein-water HBs network in the above conditions, we used Eqs 5b, 8 and 11 to calculate the time correlation C PW (t) for the intermittent HBs and the probability functions of broken hydrogen bonds N PW (t) and F PW (t), which are shown in Fig 10. The corresponding analytic expressions (four, three and two exponentials respectively) fit almost perfectly the simulation data and result to the weighted average relaxation times τ C PW , τ N PW and τ F PW , which are listed in Table 3. As it is apparent from Fig 10 and Table 3, the field slows down the kinetic functions, with 3×10 8 Hz and 3×10 9 Hz showing the slowest relaxations in C(t) and N(t), while F(t) is slower than C(t) and N(t) at all frequencies, except f 2 . We recall that these frequencies correspond to the characteristic relaxations attributed to the hydration water and mobile protein segments. The clearest improvement in fit for reaction flux k(t) is observed at 3×10 8 Hz after taking into account the F correlation, as shown in Fig 11. In fact, it is not easy to interpret the differences in the populations of HBs involved in C(t), N(t) and F(t) dynamics in relation to the applied field. The most tangible number regarding HBs kinetics is the HB life time calculated as z PW = 1/k 1 which is presented in Table 2. The largest acceleration effect of the field on HB breaking is observed at f 2 . One may expect the kinetics of PW HBs to correlate with the induced by the field positional fluctuations along the protein. The frequency with the largest impact on the protein ΔRMSF is f 1 , while that with the lowest one is f 4 . However, HBs kinetics do not correlate directly with ΔRMSFs, because HBs breakings are a combined result of protein site fluctuations and water molecule rotations each with different frequency dependence. S1 water hydrogen bonds. In this section we follow the previous analysis in order to study the kinetics of HBs between the water molecules of the first hydration shell (S1). As in the PW case, after an initial breaking of HBs between a target S1 water molecule with its neighbors, some of the neighboring molecules can be found at a distance slightly greater that 3.5 Å, still having the possibility to reform directly the initial HB after few simulation steps. The time correlation functions for the intermittent hydrogen bond C WS1 (t) and the probability of broken hydrogen bonds N WS1 (t) and F WS1 (t) are depicted in Fig 12. The corresponding average life times are listed in Table 4. A first observation is that WS1 kinetics are faster than those of PW. It seems that breaking of hydrogen bonds between S1 water molecules is 5-10 times faster than between protein and water. The N(t) kinetic is about 1.5 times faster in the WS1 system, resulting to an inverse relationship between C(t) and N(t). Apparently, rotational jumps of water molecules responsible for HB breaking are more frequent in inter-water HBs than in PW, because the protein sites are less mobile. We recall that a subset of water molecules involved in PW HBs are involved also in inter-water HBs. As it is apparent from Fig 12, the values of the above correlation functions vary with frequency but not significantly. Nevertheless, the slowest C(t) and N(t) relaxations occur at 3×10 9 Hz, whereas F(t) at 3×10 8 Hz. Again, the incorporation of F(t) improves the fits. Concerning the HB life times z WS1 , there is no noticeable effect of the field (Table 2).

Conclusions
This work has demonstrated that the dielectric spectrum of lysozyme calculated using an oscillating electric field at different frequencies, is well described by the Fröhlich-Kirkwood approach in the framework of the linear response theory. This result could be used as a criterion for choosing the maximum acceptable field strengths for similar calculations. The  Table 4. Relaxation times τ of S1 HBs correlation functions C WS1 (t), N WS1 (t) and F WS1 (t) at four characteristic frequencies and at zero field as well as hydrogen bond breaking rate constants k 1 , k 2 and k' 2 . Rate constants k 1, k 2 and k' 2 are adjusted to satisfy rate Eq (10) which includes the F(t) correlation of the free water molecules.
doi:10.1371/journal.pone.0169505.t004 Table 5. Densities of water in S1 and excess densities percentage Δ(S1:Bulk) at different frequencies and zero field conditions. Points represent k(t) as -ΔC/Δt calculated from simulation data. Blue line represents the best fit that satisfies Eq (7), while the red one is a fit including the relaxation of F(t) correlation.

Density
doi:10.1371/journal.pone.0169505.g011 Table 3. Weighted average relaxation times τ PW of PW HBs correlation functions C PW (t), N PW (t) and F PW (t) at four characteristic frequencies and at zero field as well as hydrogen bond breaking rate constants k 1 , k 2 and k' 2 . observed frequency dependent non-thermal effects on structure as well as on protein hydrogen bonds, as were described in the present study, provide new insight which can be useful in the prediction of the protein's function under similar conditions. Regarding the first hydration shell, our findings suggest a field-dependent water density as well as frequency-dependent changes in the rates of hydrogen bond breaking and reforming between protein and water. Concerning the later, the L&C model has been extended to describe better the hydrogen bond kinetics realized in field frequencies within the microwave range. The results of this work can be used as a starting point for a more detailed interpretation of future experimental studies in such systems, in the time or the frequency domain (e.g., neutron scattering).
Supporting Information S1 Table. Table A. Values of the dielectric function. Table B. Weighted average relaxation times of PW HBs. Table C. Relaxation times of S1 HBs. Table D. Hydrogen bond breaking rate constants k 1 , k 2 and k' 2 for PW HBs. Table E. Hydrogen bond breaking rate constants k 1 , k 2 and k' 2 for S1 inter-water HBs.
(DOCX) S1 Fig. Figure A. Number of intra protein HBs at different frequencies. Figure B. Number of protein-water HBs at different frequencies. Figure C. Number of interwater HBs in S1 at different frequencies. (TIF)