Pulsed electric fields of high intensity and different durations are at the basis of electroporation, an extremely promising technology developing since 90’s, with important applications in the fields of: biomedical science [1-3]; biotechnology with useful implications in food preserving [4]; nanotechnology with outcomes for bio-sensing [5] and protein analysis for binding characterization [6].

It is well known that intense electric fields can profoundly affect the physical properties and reactivity of dissolved molecules. Based on experimental results, they have been proven to play an essential role in proteins behavior, including folding, molecular recognition, and catalytic functions [7, 8]. Although one of the best acknowledged physical mechanisms at the basis of ms- and µs-pulsed electric fields is a significant increase of cell membrane permeability due to pores creation [9], when the pulse duration goes down to the nanosecond scale, and thus signal intensities reach tens of MV/m, different mechanisms can be postulated: in particular a possible effect on protein denaturation and folding [10].

Given this framework, the investigation of protein-unfolding effects induced by different waveforms of high intensity electric field pulses becomes extremely interesting. In fact the signal waveforms present in literature, at the basis of the application of nanosecond Pulsed Electric Fields (nsPEFs), are generally determined on empirical rules and often quite different from each other: monopolar pulses, bipolar ones, different pulse duration, different intensities, different rise and fall times [5, 11, 12]. Thus a rigorous study of protein-induced effects would be useful for a real comprehension of the mechanisms underlying possible protein denaturation.

A powerful tool to come up with this requirement can be provided by fully atomistic molecular simulations, which furnish affordable and reliable outcomes, a sort of “virtual experiments”, comparable to the ones provided by real experiments. Recently these techniques have been suggested as the best approach for shedding light on the different hypotheses on the interaction between electromagnetic fields and biological targets [13, 14], for which a complete theoretical investigation is still lacking. In particular fully atomistic simulations have proven to be a valuable methodology to study the effect of microwaves on enzymes in solution [15, 16], to investigate the action of pulsed electric fields on biochemical reactions in confined environments with outcomes confirmed by the experimental data [17] and to provide insight in the dynamic behavior of small peptides dispersed in water solutions, under the influence of external electric fields such as those exerted by biomedical sensors and/or field-effect transistors [18].

In a recent work [19] the authors used atomistic simulations to analyze the effect of an ideal E-field rectangular pulse on the conformational properties of myoglobin, a globular protein considered in aqueous solution. They evaluated the time behavior of the dipole response observing a fast unfolding transition characterized by a very intense polarization as well as a clear effect on the amount of relevant secondary structures, both of which are considered the cause of the protein denaturation. Nevertheless rectangular E-field pulses represent a strong idealization of the perturbation applied; therefore the purpose of the present work is to comprehend if non-ideal signals can be considered as promoters of unfolding process, as well, and to investigate the role of the specific signal waveform that can be considered.

Following this idea, the fingerprints characteristics of protein unfolding are here investigated under the action of realistic signals with lower energy content with respect to the ideal rectangular; in order to comprehend how the temporal-frequency dynamics associated to these signals can affect the protein denaturation process. The signals used are: a Gaussian signals 8-ns long, resembling a monophasic pulse, and its derivative representing a Bipolar pulse waveform of about 2-ns duration, with the same intensity but with even lower energy content.

Furthermore since pulsed signals of varying shape and duration may have severe differences in their spectral content it becomes interesting to study the response of the protein system to the different frequency components of the applied signals, using the typical techniques at the basis of the time-frequency domain analysis. Such techniques, borrowed from other scientific fields dealing with strongly not stationary signals, such as acoustics [20], geophysics [21], and biomedicine [22] are here applied for the first time to the outputs of molecular simulations. Specifically, we used a spectrogram analysis of the dipolar polarization response, to evaluate how its spectral density varies with time, depending on the different pulses applied.

Comparing the effects of such signals on protein denaturation has a twofold importance. From one hand, it concurs to elucidate the specific mechanisms of interaction between proteins and the peculiar features of the spectral content of the applied electric field; on the other hand, it could help in the choice of the signals that optimize the E-field based technology effectiveness.

Simulation details

MD simulations of myoglobin in water using the GROMACS [23] package have been carried out. The simulated system consisted of a 6-nm side rectangular box, where a single myoglobin and 6740 Single Point Charge (SPC) [24] water molecules resulting in a typical density of 1000 kg/m³. In order to properly describe myoglobin physiological behavior, it was necessary to simulate a box of water molecules large enough to reproduce both the first hydration shells and a significant amount of bulk water. The myoglobin is made of 153 residues and the Heme group (the overall charge of the system is zero).

Following an energy minimization and subsequent solvent relaxation, the system was gradually heated from 50°K to 300°K using short (typically 60 ps) MD simulations. As reference trajectory, i.e. the trajectory from which the physiological properties of myoglobin were evaluated and, especially, whereby data obtained in presence of the electric field perturbations will be compared, the one obtained in [19] has been used.

Possible effects induced by a high exogenous electric field have been evaluated by introducing two different E field signals, a Gaussian pulse (time length=8 ns) and a Bipolar pulse (time length=2 ns), whose spectrograms are reported in the upper panels of (Figure 2), Each signal has been considered with two different intensities of 10⁹ and 10⁸ V/m, acting in the simulation box as explained in [19]. More precisely the application of the electric field takes place in continuity at the last frame of the unexposed simulation, thus allowing a direct evaluation of the characteristic response over time of the system due to the switch on of the exogenous perturbation.

In particular, the introduction of the analytical form of the Gaussian and the Bipolar signals inside MD simulations has been pursued by modifying the source code of sim_util.c in the GROMACS package [23]. It is worth noting that the tinfoil Boundary Conditions have been adopted, i.e. considering negligible the effect of reaction fields on the actual simulation box and hence considering the effective electric fields acting within the simulation box as the ‘real’ electric field externally imposed by the operator.

All the trajectories were performed in the NVT (Number of particles, Volume and Temperature are constant) ensemble using an integration step of 2 fs and removing the myoglobin centre of mass translation but with no constraints on its related rotation. The temperature was kept constant at 300°K by the Berendsen thermostat [25] with the relaxation (τ) equal to the simulation time step, hence virtually equivalent to the isothermal coupling [26] which provides consistent statistical mechanical behavior. All bond lengths were constrained using LINCS algorithm [27]. Long range electrostatics was computed by the Particle Mesh Ewald method [28] with 34 wave vectors in each dimension and a 4^th order cubic interpolation. The ffG43a1 force field [29] parameters were adopted.

Numerical observables

A measure of the effects produced by the applied field is the induced protein polarization response. Field effects on protein internal behavior are more interesting rather than the trivial re-orientation of the whole protein (i.e. rotational effects); therefore the protein dipole moment has been expressed in the internal Cartesian reference of frame as explained in details in [17, 19]. Due to the non-stationarity of the applied signals it seemed worth to analyze how the frequency content of the protein dipole moment varies in time following the pulse signals. Therefore, beside the standard analyses in the frequency domain, a time-frequency processing has been performed using the spectrogram technique. The spectrogram is one of the most common time-frequency representations, able to describe the spectral content of a signal over time. The spectrogram is obtained as the square modulus of the short time Fourier transform algorithm. This algorithm implies the windowing of the temporal signal under analysis and the fast Fourier transform of each time sequence. Here, temporal sequences 5.12 10-11 s long, with 50% overlap between adjacent segments has been chosen. This allowed us to obtain a frequency resolution of 1.95 GHz and a time resolution of 0.256 ns, useful to follow the changes in time of the frequency content of the traces during the exposure. Moreover, in order to assess possible unfolding transitions induced by the external signals considered, we analyzed both the time course of hydrogen bonds disruption and/or formation (inside the protein and between the protein and the solvent) and the secondary structure time-evolution along the MD trajectories [30].

In Figure 1 the spectrogram of the protein dipolar response relative to the reference unexposed simulation is presented, showing that the power spectral density of the protein vibrational modes is mainly below 5 GHz. This indicates that such vibrational modes have a quite broad spectral content in physiological conditions, between 2 and 5 GHz. Variations in time of the frequency content evident from (Figure 1) are related to the intrinsic randomness of the protein dipolar response. An observation time of 20 ns has been considered both for the reference unexposed simulation and for the simulations with the electric signals applied, hence allowing a detailed description of electric perturbation effects during signal transitions (on-off and off-on).

Figure 1: Spectrogram of the protein dipole with no field applied.

Following the idea that pulses of lower energy content with respect to the ideal square one (1/4 for Gaussian and 1/16 for the Bipolar pulses respectively) but with higher frequencies within their duration, can couple with the intrinsic spectral content of the protein dipole previously shown, in the upper panel of (Figure 2) authors compared the spectrograms of the Gaussian pulse (B) and the Bipolar pulse (C) with the one which has been used in [19] and reported here as a reference for exposure: the rectangular electric pulse (A). The first observation is that the parameters chosen for the spectrogram calculation are suitable to follow the fast evolution in time of the spectral contents of the considered signals, thus they are suitable even for the dipole moment analysis. As expected, for the rectangular pulse, the highest frequency content is in correspondence of the pulse transitions, whereas, during the plateau, the spectrum is completely confined below 2 GHz. Conversely, the power density has a frequency content up to 4 and 6 GHz for the Gaussian and the Bipolar pulses respectively, more uniformly distributed in time, even if a higher frequency content during a shorter time interval is observable for the Bipolar pulse.

Figure 2: Spectrograms of electric field signals applied are reported in the upper part of the figure. The lower part reports the spectrograms of myoglobin dipole moment for: ideal rectangular pulse (A); Gaussian pulse (B); Bipolar pulse (C) respectively.

More in particular, for the rectangular pulse, during the on-off and off-on transitions high frequency components of the dipolar response swiftly arise (5-15 GHz), showing that a fast polarization process is taking place, able to follow the frequency components of the pulse. This transition reaches a steady state within a few hundreds of picoseconds and, after the pulse is switched off, the system does not fully relax toward the physiological state but it remains significantly polarized. This effect is clearly noticeable comparing (Figure 2) (lower panel A) with (Figure 1), reporting the reference condition. Even after the pulse switch-off, the power spectral density of the dipole moment presents a residual polarization higher than the reference level. When the perturbing electric field is a 8 ns Gaussian pulse, a similar behavior is present: the signal is still able to induce a fast polarization although the vibrational modes involved are of a lower frequency during the transitions. Results are still comparable to the previous one both during the pulse and the remaining time of the simulation (signal switched-off), as expected considering the frequency content of the Gaussian signal. Finally, the effect of a zero-mean shorter E-pulse, with lower energy content but with comparable frequency components during the duration of the signal, has been evaluated with a 2 ns Bipolar pulse (Figure 2, panel C). In this case, according to the higher frequency content of the applied field, the vibrational modes solicited are in a higher frequency range with respect to the Gaussian pulse (up to 5-6 GHz); interestingly, even with such an ultra-short signal the polarization process seems present, although a bit less intense than in the previous cases.

In summary lower energy signals seem to be able to excite the same frequencies during the time evolution of the protein dipole as the ones coupled to the ideal rectangular pulse suggesting their potential ability to induce the same effect in protein denaturation process.

To deepen the understanding of unfolding process of the protein, an insight on hydrogen-bonding of the whole system simulated has been evaluated following its evolution in time. In particular the Protein- Protein and Protein-Solvent mean number of hydrogen bonds are provided at different exposure conditions. (Table 1) clearly indicates that the three pulses drastically modify such numbers, lowering the protein-protein H-bonding and increasing the protein-solvent ones. More in particular the Gaussian pulse gives essentially the same effect of the (ideal) rectangular pulse [19], while looking at the effect of the short Bipolar pulse the effect is slightly less than the previous ones. Such an effect can be a consequence of the rearrangement of the aprotein conformation essentially due to a severe disruption of the -helix elements, which represent almost the 70% of the polypeptide, as suggested in [19].

Table 1: Protein-Protein and Protein-Solvent mean number of hydrogen bonds at different exposure conditions. When indicated the noise corresponds to the standard error.

In Figure 3 it is shown the time-course of the alpha-helix elements for the two different 109 V/m nsPEF simulations, around the instant of the onset of the electric pulse (in continuity with the reference unexposed trajectory, i.e. the one describing the physiological state of myoglobin protein) to better understand possible effects induced by the electric field just after the pulse application. Data shown provide evidences to this point, with a clear indication about myoglobin response during and after the pulse application. From the figure it is possible to appreciate the irreversible effect on protein secondary structure induced by the Gaussian pulse (Panel A of Figure 3) with essentially the same effect exerted by the ‘ideal’ rectangular pulse, although here the Gaussian nsPEF is above a value of 5*108 V/m, that can be considered a threshold for the effect, for only 2 ns. While taking into account the effect induced by the Bipolar pulse (panel B of the Figure) the relationship between pulse duration and secondary structures reorganization becomes even more pronounced: it is evident that just after the pulse crosses the threshold for polarization, the alpha-helix element starts to swiftly decrease, until the positive phase of the signal remains above 5*108 V/m, i.e. for about 300 ps. Interestingly, just when pulse intensity becomes sub-threshold, the alpha-helix elements reach an “intermediate” condition indicating that the positive peak alone is not able to induce complete unfolding. Nonetheless, when the Bipolar pulse reaches its negative peak, the alpha-helix structures undergoes a further decrease towards an irreversible unfolded state, as apparent from the secondary structures profile after the Bipolar pulse switch off (after 4.5 ns, see inset of Figure 3, panel B).

Figure 3: Red lines in the two panels correspond to the analytical waveforms of Gaussian pulse (A) and the Bipolar pulse (B); blue triangles indicate myoglobin alpha-helix elements during the electric pulses application; in the insets a magnification of the time interval when the E signal is applied and the corresponding effect on alpha-helix elements is presented, with a visual indication of E-field threshold for alpha-helix disruption (green lines).

From a comparison between the effects exerted by the Gaussian and the Bipolar nsPEFs, one can derive a sort of threshold for the E-pulse duration in order to achieve an irreversible denaturation. From the inset of Figure 3, panel A, a qualitative estimate of this time duration can be derived from the steep transition from a pre-threshold state (˜3.3 ns from the Gaussian pulse switch-on) till the minimum number of alpha-helix elements is reached (˜4.5-5 ns from the Gaussian pulse switch-on). This period is approximately 1-1.5 ns, which results in good agreement with the picture of literature for experimental nsPEFs reported in [11].

The effects produced by 108 V/m electric pulses have been investigated and, no significant effects on secondary structure properties have been obtained for all the signals, hence indicating no unfolding transitions in our simulations, similarly to what previously found [19]. This is due to the fact that these lower intensity electric fields affect myoglobin dipole moment, with a dipole transition occurring with a much slower kinetics (> 20 ns, [19]) compared to the corresponding transition occurring when a 109 V/m nsPEF is used. Thus, given the time durations of the signal adopted, no possible signal transduction can be appreciated from the corresponding spectrograms.

Molecular dynamics simulations represent a useful tool to investigate the interaction of EM fields with protein molecules. The study of polarization processes is able to provide information on major effects as denaturation or large conformational changes, in order to clarify how myoglobin activity can be controlled. In this paper we used solvated Mb to investigate in details structural effects due to Gaussian or Bipolar electric pulses with very short durations. The perturbation of the physiological state of the protein via nsPEFs produce a dramatic structural rearrangement within few hundreds of picoseconds from the pulse application, i.e. with a kinetics significantly faster, when compared to experimental data on the same protein, obtained with “standard” perturbing factors such as temperature, pressure or pH changes [32-34]. In the latter cases, in fact, myoglobin unfolding rates range from milliseconds up to seconds (or even minutes), depending on the intensity of the perturbation ; conversely in our case myoglobin undergoes a significant and fast unfolding transition only when the applied perturbing electric signal exceeds a certain amplitude.

Data have shown that the rearrangements of alpha-helix elements are sensitive to the application of either a monophasic pulse (the Gaussian one) or a biphasic one (the Bipolar pulse). The degree of unfolding induced by the Gaussian pulse is perfectly in line with the one, already published by the authors, in [19]. With reference to the Bipolar pulse, after producing a sort of intermediate state, it induces a final plateau lower than the one produced by the monophasic pulse. The interpretation of such a result seems more related to the lower energy content transferred from the signal to the protein and to the consequent elicitation of intermediate states of unfolding, rather than to a possible reversible process. In fact quite in agreement with the literature [34], which predicts slow refolding rates yet in the millisecond-second time range, and in line with results already shown in [19], the signal switching off does not induce a recovery of the protein in time durations of the order of tens of ns. Therefore it may be assumed that the fast transition between the positive and negative intensity of the Bipolar pulse shifts the protein away from the intermediate value up to the final one. In order to follow a possible reversible process after the signal switching off, longer observation times are needed, actually not easily handled by full atomistic molecular dynamics simulations.

In both cases regardless of the pattern of the signal applied, provided that it is alive for few hundreds of ps and it has a frequency content of some GHz, the protein is getting polarized and deformed following the forcing signal and, interestingly, it remains such even when the signal is removed.

These outcomes may be of particular interest since they permit a deeper comprehension of the mechanisms related to the promotion of unfolding process, clarifying the role of the specific signal waveform used, as well as the development of challenging nanotechnology applications.