Theoretical–computational modelling of the electric field effects on protein unfolding thermodynamics

Abstract

In this paper we present a general theoretical–computational approach to model the protein unfolding thermodynamics response to intense electric fields. The method proposed, based on atomistic simulations requiring limited computational effort, seems very promising to predict the unfolding thermodynamics field dependence, shedding light on the mechanisms involved. Application to myoglobin indicates a well defined field interval for a significant unfolding–refolding equilibrium with a melting field intensity ranging from 5.5 × 10⁷ to 6.0 × 10⁷ V m⁻¹ according to the protein–solvent system geometrical shape, suggesting a similar behaviour for other globular proteins.

---

Introduction

Despite the huge amount of literature on protein folding–unfolding transitions^1–7 the inherent mechanisms and their relations to the different observables experimentally monitored remain elusive and hence still a challenge for the scientific community. Interestingly in recent years, beyond the traditional studies on protein unfolding and refolding based on temperature,^8–11 pH variations^12,13 and denaturants,^14–18 the use of exotic unfolding conditions, such as pressure increase^19–23 and, very recently, intense electric fields,^24–26 has provided new stimulating perspectives to understand and possibly utilize protein unfolding processes. In particular, the effect of electric fields on protein unfolding thermodynamics has not been yet theoretically addressed and only, at the best of our knowledge, very few experimental works exist in literature on proteins response to electric fields. This is quite surprising, considering the recent technological developments based on the use of intense ultra-short electric pulses^27–30 and the increasing interest on electroporation techniques for nanomedicine, which have shown that nanosecond pulse electric field (nsPEF),^31–34 with the intensity of the order of MV m⁻¹, are able to modulate intracellular structures and functions with, moreover, a direct interaction with the genetic material.³⁵ The study of the interaction mechanisms between nsPEFs and the biological targets, due to the nanoscopic time-spatial resolution involved, can definitely take advantage of a theoretical characterization at molecular level.^36–42

In recent works^24,25 in vitro experiments showed that nsPEFs with 10⁶ ÷ 10⁷ V m⁻¹ intensities can have direct effects on enzyme activity, although circular dichroism showed intactness of the secondary structures of the enzyme²⁵ thus indicating no unfolding transitions up to these field intensities.

Even more uncommon are those experimental works explicitly focusing on the possible protein unfolding processes induced by an external electric field. In a recent paper²⁶ the authors present a single molecule method to obtain indirect evidence of protein unfolding and infer protein unfolding processes for electric fields in the 10⁶ ÷ 10⁷ V m⁻¹ range, in contrast with the findings of nsPEFs data.²⁵

Due to the intrinsic importance of understanding the basic mechanisms involved in the protein unfolding induced by electric fields and the possible relevance of their application in nanotechnology as well as in nanomedicine, we address in this paper, following the basic preliminary results we obtained in a recent article,⁴³ the problem of constructing a robust theoretical–computational quantitative approach to model the protein unfolding thermodynamics as a function of the electric field in water–protein solutions.

Theory

General relations

The thermodynamic link between the free energy of a dielectric system and an applied external (homogeneous) electric field E₀ is given bywhere M is the system thermodynamic dipole (i.e. the mean system dipole) along the external field direction and G and A are the Gibbs and Helmholtz free energy, respectively. The previous equation can be explicitly obtained from statistical mechanics once using in the canonical or isothermal–isobaric partition function the Hamiltonian⁴⁴with [scr U, script letter U]⁰_n the unperturbed energy of the n-th quantum state or phase space position (i.e. at null external field), [scr M, script letter M]′_n the corresponding dipole along the external field direction providing the linear Stark effect or first order perturbation and [scr A, script letter A], customarily considered a constant independent of the quantum state or phase space position, providing the quadratic Stark effect or second order perturbation^44–47 (i.e. it provides the energy term due to the slight electronic density polarization induced by the applied field).

From electromagnetic theory for an ellipsoidal system with volume V reacting isotropically to the field as for typical liquid state systems, when considering the external field aligned along one of the ellipsoidal axis we have⁴⁴

where ε₀ is the vacuum permittivity, χ is the (scalar) electrical susceptibility which is independent of the system shape^48,49 and f_d is a coefficient with values between zero and one, involved in the depolarizing field, determined by the shape of the system.^50,51

Note that for such systems, homogeneously polarized in the direction of the actual electric field E (i.e. the mean electric field measured inside the system) which is also homogeneous and parallel to E₀, we have from the completely general relation between the actual field and the polarization P = M/V

P = ε₀χE (4)

the following equation,as it follows comparing eqn (3) and (4). Moreover, considering that within atomistic molecular dynamics (MD) simulations (the computational procedure we use) no depolarizing field is present, i.e. the charge density at the boundary of the infinite periodic replicas is removed, MD simulations with the presence of an applied homogeneous field necessarily correspond to needle-like macroscopic ellipsoidal systems with the major axis oriented along the field as in such a case f_d = 0 and hence E = E₀.

We will therefore consider this needle-like condition as our reference geometrical shape, with reference thermodynamic dipole given by

⁰M(V, E₀) = ε₀VχE₀ (6)

Combining eqn (1) and (3) we readily obtain along the isochore

with the integral evaluated at fixed volume. Note that we removed the f_d dependence of the Helmholtz free energy at E₀ = 0 as at null field any thermodynamic property, except possibly field derivatives, is independent of the (macroscopic) system geometrical shape.

In order to evaluate the chemical potential change due to the external field we can use eqn (1) to obtain

with n the amount of molecules of the chemical species considered and μ the corresponding chemical potential (for sake of simplicity we omit in the partial derivative subscripts the amounts of the other chemical species). Note that since our theory deals with the thermodynamic response of dielectric systems to the applied electric field, we should in principle consider only macroscopic systems made of neutral chemical species. Therefore, in order to apply the same approach to solutions involving ionic solutes, we necessarily need that the ion translational motion associated to the electric current induced by the external field be much slower than the thermal atomic-molecular motions, thus allowing neglecting the electric current effects in the model (i.e. allowing to consider the ions-solvent system as statistical mechanically equivalent to the equilibrium ensemble of an identical solution except for the removal from the Hamiltonian of the terms providing the interaction between the electric potential due to the external field and the overall molecular charges).

From eqn (8) and using eqn (3) and (5), it follows that along the isochore we have

with E(E₀) given by eqn (5).

Moreover, from eqn (7) we can also obtain the solute partial molecular Helmholtz free energy a along the isochore by

with the partial molecular volume.

From eqn (9) and (10) it then follows

where p is the pressure along the isochore and p⁰ = p(V, 0), v⁰ = v(V, 0) are the pressure and partial molecular volume at null field. Moreover, considering that χ is a geometry independent intensive property fully defined by eqn (4) and hence a function of the actual field E rather than the external field E₀, we havewith from eqn (5)

From the last equations we then obtain

readily providing at E = 0where the zero superscript indicates that the property has been obtained at null field. Note that when f_d = 0 (the reference geometry) eqn (16) and (17) provide the relations of eqn (18) and (19) at whatever field, as expected as for such a geometry E = E₀.

The weak field approximation at solute infinite dilution

In order to study a water–protein solution at the typical experimental conditions, it is convenient to consider the general and exact eqn (11) in the limit of solute infinite dilution and in the weak field conditions,⁴⁴ i.e. the thermodynamic dipole can be considered linear in the field and so in eqn (11) the susceptibility and its derivatives in n, as well as the partial molecular volume v, can be considered constant along the isochore (note that at solute infinite dilution the solvent partial molecular volume is always virtually constant along the isochore). For such conditions, with the use of eqn (18) and (19), eqn (11) becomes:with χ⁰_w and Δp_w the pure water susceptibility at null field and pressure change along the isochore, respectively. Realizing that within the weak field conditions the isochorically field independent partial molecular volume can be considered, in particular for solutes at infinite dilution, also approximately constant along the p⁰ isobar and hence μ(p⁰, E, f_d) − μ(p⁰, 0) ≅ μ(V, E, f_d) − μ(V, 0) − v⁰Δp_w (i.e. we disregard the effect of the partial molecular volume variation along the p⁰ isobar on the corresponding chemical potential change), from eqn (20) we readily obtain the (weak-field) chemical potential change along the p⁰ isobar

It is worth to note that although at solute infinite dilution the derivative of any intensive property in n at fixed pressure must vanish, its product with an extensive property in the same dilution limit, when n refers to the solute molecules, has an undefined limit value, which in general might be significantly different from zero. In order to simplify the derivations it is often possible to assume, at least as a reasonable approximation, such a product limit value still vanishing, as we did in previous papers for a different thermodynamic context.^52,53 However, in this paper where we deal with V(∂χ⁰/∂n)_p,T,E such an assumption is in general unreliable and therefore we cannot use it (in the ESI†-Appendix A† of this paper we explicitly discuss the implications and limitations of assuming a vanishing product limit value).

In order to proceed further, it is convenient to consider eqn (20) for a system in the reference geometry (i.e. f_d = 0), with n corresponding to the amount of solute molecules, when using

In fact, by inserting eqn (22) into eqn (20), with f_d = 0 and considering again solute infinite dilution, we obtain

providingwhere Δ⁰p_w is the pressure change in the reference geometry.

Moreover, considering that χ_w(ρ_w, T, E) and hence

with ρ_w the water molecular density, we can write eqn (24) asprovidingwhere the geometry independent coefficientin principle a function of temperature and density, can be considered as virtually constant in density within the limited liquid water density range, as confirmed by MD simulations of the SPC water model⁵⁴ at 55 and 60 mol l⁻¹ providing pressure variations in the weak field range (E ≤ 10⁸ V m⁻¹) basically indistinguishable within the noise (data not shown). Eqn (28) then provides, at a given temperature, the differential equation for χ⁰_w(ρ_w, T)with solutionwhere ρ⁰_w is an arbitrary reference density within the liquid range.

Using again eqn (22) into eqn (20) with now f_d ≠ 0 we obtain the generalization of eqn (24) for whatever ellipsoidal shape

readily providing, when expressing Δp_w = Δ⁰p_w + (Δp_w − Δ⁰p_w),

Hence, from eqn (24) and (27) we have

providing the pressure change along the isochore for an arbitrary ellipsoidal shape. It is worth noting that eqn (33) shows that the presence of the depolarizing field increases the pressure change and for a system with a large susceptibility as liquid water, such an increase can be very large as f_d → 1.

Using eqn (21) and (33) we then have

providing the solute chemical potential change as a function of the actual field along the isobar, withcorresponding to the (null field) susceptibility change due to inserting a single solute molecule, at fixed volume, in a pure water (macroscopic) system. Interestingly, comparison of eqn (34) with eqn (20) clearly shows that for a given actual field the geometrical shape effects on the chemical potential are present only in isobaric conditions. From eqn (20), (34) and (35) we then obtain the chemical potential change due to the unfolding process μ_D − μ_Nwhere χ⁰_D, χ⁰_N are the (null field) susceptibilities when inserting either the unfolded (denatured) or the folded (native) protein molecule, v⁰_D and v⁰_N are the corresponding protein partial molecular volumes and from the definition of p⁰ we have μ_D(p⁰, 0) − μ_N(p⁰, 0) = μ_D(V, 0) − μ_N(V, 0).

Finally, it is worth to consider that although in eqn (21), (34) and (36) we disregarded the partial molecular volume variation along the p⁰ isobar as its effect on the chemical potential change should be negligible (in particular for the highly diluted solutes), the solvent molecular density along the p⁰ isobar cannot be in general virtually constant, as occurring along the isochore, as a significant pressure change between the isochore and the isobar can be present (see eqn (33)). In fact, the pressure of a pure liquid water system p_w, identical to the pressure of aqueous solutions with solutes at infinite dilution, at each actual field value for a given ellipsoidal shape can be expressed by

p_w(ρ_w, E, f_d) = p_w(ρ_w, 0) + Δp_w(ρ_w, E, f_d) (38)

thus providing for the state points along the p⁰ isobar

p⁰ = p_w(ρ_w, 0) + Δp_w(ρ_w, E, f_d) (39)

Therefore, considering that within the water liquid range the isothermal compressibility

is roughly constant hence implying thatwhere ρ_w^(p0) is defined by p⁰ = p_w(ρ_w^(p0), 0) and we used the fact that liquid water is highly incompressible and thus , by means of eqn (30), (33) and (41) we can solve eqn (39) explicitly obtaining the liquid water density corresponding, along the p⁰ isobar, at each E value (within the weak field range) for a given ellipsoidal shape.

Computational methods

We carried out NVT MD simulations of pure SPC⁵⁴ water model systems at 300 K for different density conditions, ranging from 55 mol l⁻¹ to 60 mol l⁻¹, in order to characterize the SPC susceptibility and compare such data with the experimental water electric susceptibility data available in literature.⁵⁵ According to the general theory described in a previous paper,⁴⁴ the electric susceptibility of pure, liquid state, water model systems can be obtained by using a uniform distribution model for the dipole fluctuations, accurately reproducing the field dependence of the system thermodynamic dipole (i.e. the mean dipole) as obtained by MD simulations⁴⁴ (note that MD systems correspond to needle-like ellipsoidal systems with the applied field along the major axis, our reference condition, as no depolarizing field is present in the simulations and hence E = E₀). Such a model properly reproduced the SPC dipole-field behaviour in the whole density range we considered and, interestingly, the null field susceptibilities obtained are virtually indistinguishable from the corresponding values provided by using the weak field approximation (i.e. assuming a linear behaviour of the thermodynamic dipole equivalent to using the Gaussian distribution model⁴⁴) up to E = 10⁸ V m⁻¹, indeed confirming our previous results⁴⁴ showing that the weak field approximation is highly accurate for E ≤ 10⁸ V m⁻¹ (relative deviations between the two sets of susceptibilities were always within 1%).

The pure SPC system at 55.32 mol l⁻¹ and 300 K, corresponding to the typical experimental liquid water condition, also served as reference system to insert a single protein molecule at fixed volume. In this case the simulated system consisted of a cubic box (about 6 nm side) where we placed a single sperm whale myoglobin and 7202 SPC water molecules.

The myoglobin is made of 153 residues plus the heme group (the overall protein charge is zero with hence no need of counter-ions in the simulations). 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. A first trajectory was propagated up to 50 ns in the NVT ensemble using an integration step of 2 fs. For both the pure SPC and SPC–myoglobin solution systems the temperature was kept constant at 300 K by the isothermal coupling⁵⁶ which provides a consistent statistical mechanical behaviour. All bond lengths were constrained using the LINCS algorithm.⁵⁷ Long range electrostatics were computed by the Particle Mesh Ewald method⁵⁸ with 34 wave vectors in each dimension and a 4th order cubic interpolation. The ffG43a1 force field⁵⁹ parameters were adopted. Short range interactions were evaluated within a 1.1 nm cut off radius.

Once obtained an extended equilibrated-unexposed trajectory we obtained a set of exposed trajectories by following the same approach used in a recent paper,⁴³ i.e. by means of simulations with a stationary exogenous electric field of increasing intensity (from 10² to 10⁹ V m⁻¹). In agreement with our previous results,⁴³ we obtained a clear unfolding transition followed by a stable unfolded state of the protein at E ≥ 5.0 × 10⁸ V m⁻¹. Starting from the equilibrated folded and unfolded myoglobin (the latter obtained by the MD simulation at the highest field intensity) we also performed a set of NVT MD simulations at the same SPC density with null, 5.0 × 10⁶, 10⁷, 5.0 × 10⁷ and 10⁸ V m⁻¹ stationary exogenous field for the folded conformation and with null, 10⁸, 5.0 × 10⁸, 10⁹ V m⁻¹ for the unfolded conformation, now using a large simulation box (9 nm side and 24291 SPC molecules) ensuring a minimum distance between protein atoms and the simulation box faces always significantly larger than 1 nm, thus including both the protein hydration shell⁶⁰ and a significant amount of bulk solvent molecules according to the ≈1 nm thickness of the hydration shell reported in literature.⁶¹ It is worth noting that within the production runs of these latter simulations (each with a 10–20 ns time-length) the unfolded protein remained stable and no refolding transitions occurred (i.e. the refolding kinetics is too slow). Finally, two further simulations starting from the equilibrated folded and unfolded myoglobin at null field with the same large box volume but with a slightly modified number of SPC molecules, were performed in order to have the SPC–myoglobin system at the same pressure of the pure SPC reference system. For all the SPC–myoglobin simulations within the large simulation box all the simulation parameters (except obviously the SPC amount and the box size) were identical to the ones described for the other simulations.

In Fig. 1 it is reported the thermodynamic dipole for the myoglobin solution when the protein is either folded or unfolded, at different exposure conditions (we always evaluated the simulation box mean dipole considering the protein in the center of the box, i.e. re-centering the protein at each MD frame, in order to prevent any protein molecule disruption due to the periodic boundary conditions). In particular for the folded case, only data up to 10⁸ V m⁻¹ are present, since at higher field intensities unfolding transitions occurred within the simulation time-length. Conversely, for the unfolded state, at each electric field condition the unfolded conformation was stable and, therefore, we were in principle able to collect mean dipole moment data over the whole field intensity range considered. However, the unfolded myoglobin solution simulations with field intensities lower than 10⁸ V m⁻¹ (i.e. field intensities unable to keep the protein dipole closely aligned along the field direction) might be unreliable due to the slow rotational kinetics of the protein with thus a possible relevant overestimation of the system mean dipole (note that the myoglobin unfolded conformation is characterized by a large intrinsic dipole). Moreover, the myoglobin unfolded conformation as provided by the 10⁹ V m⁻¹ MD simulation, although reasonably representative of the equilibrium unfolded state for E ≥ 10⁸ V m⁻¹ (at 5.0 × 10⁸ V m⁻¹ a fast unfolding transition is apparent within the simulation time-length), could be not the proper conformation to be used for the unfolded state simulations at lower field intensities, in particular close to the null field condition (note that for both the folded and unfolded myoglobin solutions we always considered an exactly zero mean dipole moment at null field, thus removing the sampling noise). Therefore, due to the possible relevant dipole noise and conformational uncertainty in the protein unfolded state simulations with low field intensities, we did not consider in our calculations the unfolded myoglobin MD simulations with E < 10⁸ V m⁻¹.

Fig. 1 Mean dipole moment along the field direction as a function of the electric field intensity as obtained by SPC–myoglobin simulations. The inset magnifies the weak field range.

In order to evaluate the susceptibility change between the unfolded and folded myoglobin solutions to be used in eqn (36) and (37) for obtaining μ_D − μ_N, we made use of the weak field approximation in the range 0 ÷ 10⁸ V m⁻¹ where both the unfolded and the folded conformations are stable within the simulation time-length (we assumed that within such a field range either for the folded or the unfolded state myoglobin no conformational/structural transitions possibly occurring can significantly change the mean dipole linear behaviour of the corresponding solution). Note that for the unfolded state condition, the use of the uniform distribution model⁴⁴ to evaluate over the whole field range the system susceptibility provided a value almost identical to the susceptibility as obtained by the weak field approximation (relative deviation within 3%).

Moreover, we expect that the term ε₀V(χ⁰_D − χ⁰_N) as obtained by the weak field approximation be even closer to the corresponding exact value as the slight errors of the unfolded and folded susceptibilities as obtained by the weak field approximation, are likely to cancel out in the difference.

Finally, it is worth to remark that although the use of nanoscopic simulation boxes to estimate thermodynamic derivatives in the solute molecular number, at fixed volume, could be affected by a systematic error due to the possible relevant pressure increase (especially when dealing with large solutes), in the present case where we only need to estimate the folded to unfolded state susceptibility change we can rather safely assume such nanoscopic effects to be negligible, at least when using a simulation box including the protein hydration shell and a significant amount of bulk solvent molecules (see ESI-Appendix B†). In fact, we compared the myoglobin structural behaviour as provided by the null field MD simulations at fixed SPC reference density (i.e. corresponding to inserting the protein into the SPC box at fixed volume) with the corresponding structural distributions as obtained by the simulations of the same SPC–myoglobin system at null field with the same box volume but with a slightly modified number of SPC molecules in order to have the same pressure of the pure SPC reference system and hence mimicking the insertion of the protein into the SPC box at fixed pressure. The two sets of simulations provided Gaussian-like equilibrium distributions of the protein structural sub-states corresponding to different protein excluded volumes, for both the folded and unfolded conformations, almost identical (distribution parameters indistinguishable within the noise and/or with a relative deviation within 1%) and, moreover, the protein mean excluded volume change for the folded to unfolded state transition indistinguishable, within the noise (see Table 1). Such results indicate that protein chemical potential variations due to the structural transitions, in particular for the folded to unfolded state transition, were essentially unaffected by the nanoscopic size effects and hence should be considered as thermodynamically converged within a reasonably small error.

Table 1 Protein excluded volume distribution parameters for the folded and unfolded conformations, as obtained from the null field simulations at fixed reference SPC density or pressure. In addition, we also show the mean excluded volume change for the folded to unfolded state transition for both the simulation sets. The noise corresponds to two standard errors, evaluated by using 5 different portions of the total trajectory

SPC–myoglobin system	Excluded volume for the folded conformation		Excluded volume for the unfolded conformation		Folded to unfolded excluded volume variation
	Mean value (nm^3)	Variance (nm^6)	Mean value (nm^3)	Variance (nm^6)	Mean difference (nm^3)
Reference SPC density	26.8 ± 0.06	0.89 ± 0.04	27.65 ± 0.09	1.00 ± 0.06	0.85 ± 0.11
Reference SPC pressure	27.09 ± 0.10	0.97 ± 0.04	27.80 ± 0.06	0.98 ± 0.03	0.71 ± 0.12

---

Results

In order to test our theoretical model we first considered the susceptibility of pure liquid water systems, comparing the available experimental⁵⁵ and simulation data (for the latter we used our set of SPC simulations at different densities) with the expected linear relation between the null field water susceptibility and density predicted by our theory (i.e. eqn (30)).

From Fig. 2and 3 it is evident that, within the liquid range, both the experimental water data as well as the SPC simulation data are fully consistent with the predicted linear behaviour provided by eqn (30) of our theoretical model. Interestingly, the linear regression of the experimental water data at 300 K provides slope and intercept rather close to the ones obtained by the linear regression of the corresponding simulation data of SPC at the same temperature, showing that the SPC model, within the liquid state range, reproduces reasonably well the liquid water dielectric behaviour.

Fig. 2 Experimental liquid water susceptibility at different temperatures as a function of the density ratio ρ_w/ρ⁰_w (circles) and corresponding linear regressions (solid lines). The reference density ρ⁰_w = 55.508 mol l⁻¹ is taken according to “static dielectric constant of water and steam” (1980).⁵⁵

Fig. 3 Pure SPC susceptibility as a function of the density ratio ρ_w/ρ⁰_w as obtained by MD simulations at 300 K (circles) and corresponding linear regression (solid line). The reference density ρ⁰_w = 55.508 mol l⁻¹ is taken according to “static dielectric constant of water and steam” (1980).⁵⁵

In Fig. 4 the distributions of the helical secondary structure content (including the alpha helix, 3-helix and 5-helix) are shown for both the folded and unfolded state myoglobin simulations at different electric field intensities. The figure clearly demonstrates that for the folded myoglobin no significant secondary structure disruption, with respect to the null-field condition, is present when the applied electric field is up to 10⁸ V m⁻¹. Therefore, such results ensure that within the whole weak-field range our folded myoglobin simulations are suitable to properly describe the folded state condition. Interestingly, the same figure also shows that for the unfolded myoglobin simulations (at 10⁸ V m⁻¹ and 10⁹ V m⁻¹) about 80% of the secondary structure content is lost, with no relevant variation due to the electric field intensity change (in the ESI,† Appendix C,† we show for the same field intensities considered in Fig. 4 the corresponding single residue secondary structure timelines,⁶² evidencing that for the folded myoglobin no relevant secondary structure changes are present).

Fig. 4 Distributions of the myoglobin helical secondary structure content (including alpha helix, 3-helix and 5-helix) for both the folded and unfolded myoglobin simulations as a function of the applied electric field.

It is worth to note that within the whole weak-field range no relevant conformational–structural transitions can be detected for the folded myoglobin, as illustrated by Fig. 5 where we show the distributions at different electric field conditions, up to 10⁸ V m⁻¹, of the heme plane orientation with respect to the myoglobin major geometrical axis⁴³ (in the ESI, Appendix C,† we also show the time course of the radius of gyration⁶³ and the solvent accessible surface area⁶³ for the same field intensities considered in Fig. 4).

Fig. 5 Distributions of the projection of the unit vector orthogonal to the heme plane over the myoglobin major geometrical axis, as a function of the electric field, for the folded myoglobin simulations. The protein geometrical axes were obtained by means of diagonalisation of the 3 × 3 geometrical covariance matrix.⁴³

By using the susceptibility χ⁰_w(ρ_w, T), at 55.32 mol l⁻¹ and 300 K, as provided by the linear regression of experimental water data, and the value of ε₀V(χ⁰_D − χ⁰_N) as obtained by the SPC–myoglobin simulations at fixed reference SPC density, we obtained by means of eqn (36) and (37) the chemical potential change μ_D − μ_N due to the unfolding process as a function of the actual field present in the system, as shown in Fig. 6, where we used for μ_D(p⁰, 0) − μ_N(p⁰, 0) = μ_D(V, 0) − μ_N(V, 0) and v⁰_D − v⁰_N the values taken from myoglobin experimental data (52 ± 1.6 kJ mol⁻¹ for the former and about −100 ml mol⁻¹ for the latter).^64,65

Fig. 6 Chemical potential change of unfolding as a function of the actual electric field along the isobar at three different ellipsoidal geometries (solid lines). In the figure it is also shown the chemical potential change of unfolding along the reference isochore (circles). Note that the estimated error for the melting field is about 10⁶ V m⁻¹.

From Fig. 6 it is evident the geometrical effect when comparing the chemical potential change along the isobar for three significant geometries: f_d = 0, a thin needle with the major axis along the field; f_d = 1/3, a sphere; f_d = 1, a flat disk perpendicular to the field. For the needle-like geometry (see the inset of Fig. 6) the chemical potential change reaches the zero value, corresponding to the melting field (i.e. the field providing identical fractions of the unfolded and folded populations), when the actual field is about 5.5 × 10⁷ V m⁻¹, while for the sphere and even more for the flat disk a significant increase of the field is necessary to reach the same equilibrium condition (about 5.7 × 10⁷ V m⁻¹ for the sphere and about 6.0 × 10⁷ V m⁻¹ for the disk). For comparison, in the figure it is also shown the chemical potential change as a function of the field as obtained for the reference isochore, virtually coinciding with the isobaric needle-like chemical potential change. It is worth to note that the estimated relative error of the field variation of the chemical potential change of unfolding close to the melting field is about 2% which, combined with the experimental noise of the null field chemical potential change of unfolding, provides an estimated total error for the melting field of about 10⁶ V m⁻¹.

In Fig. 7 it is shown the unfolded state equilibrium fraction for the isobaric folding–unfolding equilibrium at the three geometries considered as a function of the actual electric field.

Fig. 7 Unfolded state fraction as a function of the actual electric field along the isobar at three different ellipsoidal geometries.

From this figure it results that below 4.7 × 10⁷ V m⁻¹ the unfolded state population can be neglected and beyond 6.7 × 10⁷ V m⁻¹ only the unfolded state is virtually present. In Fig. 8 we show the predicted pressure variations along the isochore, according to eqn (33), for pure liquid water at 55.32 mol l⁻¹ and 300 K (identical to the pressure variations of the corresponding water–protein system at infinite protein dilution) as a function of the actual field intensity for the three geometries considered.

Fig. 8 Isochoric pressure change of liquid water at 55.32 mol l⁻¹ and 300 K, according to eqn (33), as a function of the actual electric field for the three different ellipsoidal geometries.

The figure clearly shows that a large pressure change is induced by the geometrical transition from the needle-like to the disk shape.

Finally in Fig. 9 we show, for the same three geometries considered, the liquid water molecular density change as a function of the actual field along the isobar, indicating that within the weak field range even a large pressure change between the isochore and the isobar (see Fig. 8) results in a rather limited liquid water density variation.

Fig. 9 Liquid water molecular density versus the actual field along the isobar, for the three different ellipsoidal geometries.

Conclusions

In this paper we presented a general and robust theoretical–computational approach to model protein unfolding thermodynamics under the effects of an external electric field within the so called weak field range, i.e. providing a linear field dependence of the system mean dipole and corresponding to an actual field inside the aqueous solution up to ≈10⁸ V m⁻¹. Comparison of the prediction of our theoretical model for the liquid water susceptibility density dependence with the available experimental and computational data, showed that the level of theory used is able to capture the essential features of the dielectric response to the field and hence to provide a reliable basic description of the thermodynamic changes induced by the field in proteins aqueous solutions. Interestingly, our theoretical derivations indicate that in the presence of the electric field the geometrical shape of the protein–water system behaves as an additional state variable, beyond its effect on the actual field, providing relevant changes in the folding–unfolding equilibrium as it is varied. Application of our approach to myoglobin showed that no unfolded state population is present for field intensities below 4.7 × 10⁷ V m⁻¹, with a totally unfolded state population at field intensities above 6.7 × 10⁷ V m⁻¹. Such results suggest that for similar globular proteins the use of electric fields below 10⁷ V m⁻¹ should be unable to induce any significant unfolding, in line with the experimental evidences of no secondary structure disruption up to such field intensities, even in the presence of protein activity reduction.^24,25

Finally, given the fact that the approach outlined in this paper only requires a limited computational effort and the use of (typically) available experimental data, we expect it to be particularly suited for investigating the thermodynamics field dependence of several proteins in different physical–chemical conditions, possibly allowing a systematic quantitative study on the electric field effects on protein folding–unfolding behaviour.

References

1. J. Rösgen and H. J. Hinz, Phys. Chem. Chem. Phys., 1999, 1, 2327.
2. I. Burgos, S. A. Dassie and G. D. Fidelio, J. Phys. Chem. B, 2008, 112, 14325.
3. W. A. Eaton, V. Mũnoz, S. J. Hagen, G. S. Jas, L. J. Lapidus, E. R. Henry and J. Hofrichter, Annu. Rev. Biophys. Biomol. Struct., 2000, 29, 327.
4. B. Linder and R. A. Kromhout, J. Phys. Chem. B, 2001, 105, 6387.
5. T. Schindler, M. Herrler, M. A. Marahiel and F. X. Schmid, Nat. Struct. Biol., 1995, 2, 663.
6. F. N. Zaidi, U. Nath and J. B. Udgaonkar, Nat. Struct. Biol., 1997, 4, 1016.
7. H. J. Dyson and P. E. Wright, Nat. Struct. Biol., 1998, 5, 499.
8. R. L. Baldwin, Proc. Natl. Acad. Sci. U. S. A., 1986, 83, 8069.
9. R. Day, J. Brian, N. J. Bennion, S. Ham and V. Daggett, J. Mol. Biol., 2002, 322, 189.
10. Y. Moriyama and K. Takeda, J. Phys. Chem. B, 2010, 114, 2430.
11. G. Bellavia, G. Cottone, S. Giuffrida, A. Cupane and L. Cordone, J. Phys. Chem. B, 2009, 113, 11543.
12. D. Stigter, D. O. V. Alonso and K. A. Dill, Proc. Natl. Acad. Sci. U. S. A., 1991, 88, 4176.
13. O. O. Sogbein, D. A. Simmons and L. Konermann, J. Am. Soc. Mass Spectrom., 2000, 11, 312.
14. F. Ahmad and C. C. Bigelow, J. Biol. Chem., 1982, 257, 12935.
15. T. T. Herskovits and H. Jaillet, Science, 1969, 163, 282.
16. C. N. Pace and K. E. Vanderburg, Biochemistry, 1979, 18, 288.
17. S. Hashimoto, J. Fukasaka and H. Takeuchi, J. Raman Spectrosc., 2001, 32, 557.
18. L. Konermann, F. I. Rosell, A. G. Mauk and D. J. Douglas, Biochemistry, 1997, 36, 6448.
19. S. A. Hawley, Biochemistry, 1971, 10, 2436.
20. J. Wiedersich, S. Köhler, A. Skerra and J. Friedrich, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 5756.
21. G. Panick, R. Malessa, R. Winter, G. Rapp, K. G. Frye and C. A. Royer, J. Mol. Biol., 1998, 275, 389.
22. G. Hummer, S. Garde, A. E. García, M. E. Paulaitis and L. R. Pratt, Proc. Natl. Acad. Sci. U. S. A., 1998, 95, 1552.
23. L. L. Shen and J. Hermans, Biochemistry, 1972, 11, 1836.
24. S. J. Beebe, Bioelectrochemistry, 2015, 103, 52.
25. T. Yu and X. Fu, Bioelectrochemistry, 2015, 101, 42.
26. K. J. Freedman, S. R. Haq, J. B. Edel, P. Jemth and M. J. Kim, Sci. Rep., 2013, 3, 1.
27. C. Merla, S. El Amari, M. Kenaan, M. Liberti, F. Apollonio, D. Arnaud-Cormos, V. Couderc and P. Leveque, IEEE Trans. Microwave Theory Tech., 2010, 58, 4079.
28. C. Merla, A. Denzi, A. Paffi, M. Casciola, G. Dinzeo, F. Apollonio and M. Liberti, IEEE Trans. Biomed. Eng., 2012, 59, 2302.
29. M. Balucani, P. Nenzi, R. Crescenzi, P. Marracino, F. Apollonio, M. Liberti, A. Densi and C. Colizzi, Technology and design of innovative flexible electrode for biomedical applications, ECTC, Lake Buena Vista, FL, 2011,  DOI:10.1109/ectc.2011.5898682.
30. P. Nenzi, A. Denzi, K. Kholostov, R. Crescenzi, F. Apollonio, M. Liberti, P. Marracino, A. Ongaro, R. Cadossi and M. Balucani, Smart flexible planar electrodes for electrochemotherapy and biosensing, ECTC, Las Vegas, NV, 2013,  DOI:10.1109/ectc.2013.6575616.
31. B. L. Ibey, S. Xiao, K. H. Schoenbach, M. R. Murphy and A. G. Pakhomov, Bioelectromagnetics, 2009, 30, 92–99.
32. C. Yao, Y. Mi, X. Hu, C. Li, C. Sun, J. Tang and X. Wu, Experiment and mechanism research of SKOV3 cancer cell apoptosis induced by nanosecond pulsed electric field, 30th Annual International IEEE EMBS Conference Vancouver, British Columbia, Canada, 2008.
33. S. J. Beebe, Y. J. Chen, N. M. Sain, K. H. Schoenbach and S. Xiao, PLoS One, 2012, 7, e51349.
34. L. Chopinet and M. P. Rols, Bioelectromagnetics, 2015, 103, 2–6.
35. S. J. Beebe and K. H. Schoenbach, J. Biomed. Biotechnol., 2005, 2005, 297.
36. R. Reale, N. J. English, J. A. Garate, P. Marracino, M. Liberti and F. Apollonio, J. Chem. Phys., 2013, 139, 205101.
37. P. Marracino, A. Amadei, F. Apollonio, G. D Inzeo, M. Liberti, A. D. Crescenzo, A. Fontana, R. Zappacosta and M. Aschi, J. Phys. Chem. B, 2011, 115, 8102.
38. M. Avena, P. Marracino, M. Liberti, F. Apollonio and N. J. English, J. Chem. Phys., 2015, 142, 141101.
39. R. Reale, N. J. English, P. Marracino, M. Liberti and F. Apollonio, Mol. Phys., 2014, 112, 1870.
40. R. Reale, N. J. English, P. Marracino, M. Liberti and F. Apollonio, Chem. Phys. Lett., 2013, 582, 60.
41. P. Marracino, M. Liberti, G. d'Inzeo and F. Apollonio, Bioelectromagnetics, 2015, 36, 377.
42. F. Apollonio, M. Liberti, A. Amadei, M. Aschi, M. Pellegrino, M. D'Alessandro, M. D'Abramo, A. Di Nola and G. D'Inzeo, IEEE Trans. Microwave Theory Tech., 2008, 56, 2511.
43. P. Marracino, F. Apollonio, M. Liberti, G. d’Inzeo and A. Amadei, J. Phys. Chem. B, 2013, 117, 2273.
44. M. E. F. Apol, A. Amadei and A. Di Nola, J. Chem. Phys., 2002, 116, 4426.
45. A. Haug, Theoretical Solid State Physics, Pergamon, Oxford, 1972.
46. P. W. Atkins and R. S. Friedman, Molecular Quantum Mechanics, Oxford University Press, Oxford, 3rd edn, 1997.
47. P. W. Selwood, Magnetochemistry, Interscience, New York, 2nd edn, 1956.
48. E. A. Guggenheim, Thermodynamics, North-Holland, Amsterdam, 6th edn, 1977.
49. I. C. Yeh and M. L. Berkowitz, J. Chem. Phys., 1999, 111, 3155.
50. E. C. Stoner, Magnetism and Matter, Methuen, London, 1934.
51. B. I. Bleaney and B. Bleaney, Electricity and Magnetism, Oxford University Press, Oxford, 3rd edn, 1976.
52. M. D’Abramo, M. D’Alessandro and A. Amadei, J. Chem. Phys., 2004, 120, 5226.
53. M. D’Alessandro, M. D’abramo, G. Brancato, A. Di Nola and A. Amadei, J. Phys. Chem. B, 2002, 106, 11843.
54. H. J. C. Berendsen, J. P. M. Postma, W. F. V. Gunsteren and J. Hermans, Interaction Models for Water in Relation to Protein Hydration. in Intermolecular Forces, ed. B. Pullman, Reidel Publishing Company, Dordrecht, The Netherlands, 1981.
55. E. U. Franck, J. Phys. Chem. Ref. Data, 1980, 9, 1291.
56. D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids, Academic Press, London, 1990.
57. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Frajie, J. Comput. Chem., 1997, 18, 1463.
58. T. A. Darden, D. M. York and L. G. Pedersen, J. Chem. Phys., 1993, 98, 10089.
59. W. F. van Gunsteren, S. R. Billeter, A. A. Eising, P. H. Hunemberger, P. Kruger, A. E. Mark, V. R. P. Scott and I. G. Tironi, Biomolecular Simulation: The GROMOS96 Manual and User Guide, Hochschlverlag AG an der ETH, Zurich, Switzerland, 1996.
60. L. R. Murphy, N. Matubayasi, V. A. Payne and R. M. Levy, Folding Des., 1998, 3, 105.
61. S. Ebbinghaus, S. J. Kim, M. Heyden, X. Yu, U. Heugen, M. Gruebele, D. M. Leitner and M. Havenith, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 20749.
62. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
63. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
64. E. Bismuto, G. Irace, L. Servillo, A. Giovane and G. Colonna, Experientia, 1984, 40, 1400.
65. C. A. Royer, Biochim. Biophys. Acta, 2002, 1595, 201.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra15605j

---