{'fraction_non_alphanumeric': 0.03352010874970989, 'fraction_numerical': 0.013328470541427672, 'mean_word_length': 4.6652892561983474, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dissipative electro-elastic network model of protein electrostatics

We propose a dissipative electro-elastic network model (DENM) to describe the dynamics and statistics of electrostatic fluctuations at active sites of proteins. The model combines the harmonic network of residue beads with overdamped dynamics of the normal modes of the network characterized by two friction coefficients. The electrostatic component is introduced to the model through atomic charges of the protein force field. The overall effect of the electrostatic fluctuations of the network is recorded through the frequency-dependent response functions of the electrostatic potential and electric field at the active site. We also consider the dynamics of displacements of individual residues in the network and the dynamics of distances between pairs of residues. The model is tested against loss spectra of residue displacements and the electrostatic potential and electric field at the heme's iron from all-atom molecular dynamics simulations of three hydrated globular proteins.


INTRODUCTION
elastic network with the molecular charge distribution might be sufficient to describe the statistics of the potential fluctuations and their slow dynamics.
We propose here a model combining the distribution of molecular charge from the standard atomic force-fields with the dynamics and statistics of protein motions derived from an elastic network. To test the model, we employ all-atom Molecular Dynamics (MD) simulations of three hydrated heme proteins. Previous analysis of these data has emphasized the importance of nanosecond (ns) relaxation modes in the electrostatic fluctuations of the protein matrix and the protein-water interface (14,35). This time-scale is relevant to biological function since heme proteins are typical components of biology's energy chains transporting electrons on nanosecond-to-microsecond time-scales (7,36). The nanosecond time-scale fluctuations are also consistent with the dynamics of the global motions of the network of residues. An elastic network model naturally fits the physics of the problem. The present contribution is a first step in developing network models of the electrostatic response of hydrated proteins. Here, we focus on the electrostatics of the protein matrix only, leaving the water component of the overall response to future studies.

Elastic network model
A coarse-grained elastic network model (ENM) assigns a node to a collection of atoms reducing the computational burden of an all-atom normal-mode analysis. The typical coarse-graining is done on the level of individual aminoacids (1-bead model (37)). The position of the node is defined by the coordinates of the C α atom. The springs connecting the nodes represent the bonded and non-bonded interactions within an accepted cutoff distance (18)(19)(20) or by a distance-dependent force constant (26,33,38,39). The ENM diagonalizes the Hessian matrix H derived from a simplified Hookean potential suggested by Tirion (18). This potential, E i j = C(r i j − r 0,i j ) 2 /2, describes the elongation r i j = |r i − r j | between nodes i and j in the network characterized by one universal force constant C and the structural information stored in the equilibrium bead positions r 0,i (r 0,i j = |r 0,i − r 0, j |).
The Hessian H αβ i j is a 3N × 3N matrix representing the protein elastic energy as a quadratic form of the Cartesian displacements δr α i = r α i − r α 0,i of individual beads in the network Here, i and j run between 1 and N, α, β indicate the Cartesian projections, and summation over repeated Greek indices is assumed here and throughout. The Hessian is diagonalized by the unitary matrix U producing 3N eigenvalues λ m . This standard linear algebra formalism (18,19) yields the statistical correlator between the displacements of residues i and j in the network where β = 1/(k B T ) is the inverse temperature. If the network is characterized by a universal force constant and a cutoff, it yields a bell-shaped density of vibrational states. A refinement of that network by adopting a stronger coupling between covalently bound neighbors splits the density of states into two maxima, in better agreement with allatom calculations (24,40,41). This approximation is adopted in our present calculations: the spring constant was multiplied by a constant factor ε = 100 for covalent neighbors. In addition, a uniform cut-off radius of 15Å was used in all calculations.

Overdamped network dynamics
The standard mechanical ENM outlined above obviously lacks dissipative dynamics. The equations of motion are harmonic, implying oscillatory time correlation functions. In contrast, most correlation functions of hydrated proteins observed by scattering (8,10,42) and relaxation (43) techniques are exponential, corresponding to the overdamped (Debye) dynamics, or stretched-exponential (8). Several approaches can be implemented to incorporate dissipative dynamics into the mechanical network of beads. Langevin equations of motion for the ENM potential, within the general framework of the Lamm-Szabo formalism (44), have been suggested (34,48). This approach, and some early suggestions (45,46), still requires parameterization (done by fitting the rotational and translational diffusion coefficients from MD) and, in addition, doubles the size of the matrix to be diagonalized. We have adopted here a more straightforward formalism not requiring additional computational resources.
Instead of a harmonically oscillating network, we have assumed for each normal mode q m , diagonalizing the network's Hessian, an overdamped motion described by the dissipative memory kernel ζ(t) (47). The equation of motion for such overdamped dynamics is where F(t) is an external force. The distinction between Eq. (3) and the Langevin network (34,48) is worth emphasizing here. In the latter, uncoupled Langevin dissipative equations are first assigned to each bead of the network. However, it was noted that dissipative equations for the beads are likely to become coupled when the Langevin dynamics of beads are consistently derived by integrating out the fast degrees of freedom in the equations of motion (49). Given that the assignment of uncoupled Langevin equations to the individual beads is phenomenological from the onset, one can introduce similar phenomenology at a different level of the theory. Here, in the spirit of the standard normal-mode analysis, we introduce uncoupled dissipative equations for the normal modes (Eq. (3)). This is still a phenomenological assumption requiring further testing, but it reduces to the standard normal-mode analysis in the static limit.
When Laplace-Fourier transform (47) is applied to Eq. (3), the displacement response function for collective mode q m follows whereζ(ω) is the Laplace-Fourier transform of the friction kernel ζ(t). Extending this procedure to all normal modes of the network, one obtains the response function of bead displacements In this equation, the response function χ i j (ω), which is a rank-2 tensor, represents the displacement of residue i of the protein due to an oscillatory force with frequency ω applied to residue j and propagated through the network to residue i. This response function is the basis of the dissipative electro-elastic network model (DENM) proposed here. At ω = 0, Eq.  Electric fields E ik and E jk are produced by the corresponding charges q ik and q jk at the active site. The overall electric field of residue i at the active site, E 0i , is obtained by summation of all E ik produced by the charges q ik of the residue.

Electrostatic potential response function
The ENM has shown good performance for low-frequency structural/mechanical deformations of individual proteins and biomolecular assemblies (27). Our goal here is to supplement these low-frequency motions with atomic partial charges to model the dynamics and statistics of electrostatic fluctuations at active sites of proteins. Coarse-grained electrostatics of proteins has been addressed in a number of recent papers (50)(51)(52), mostly dealing with long-range protein assembly and interactions (53). These approaches also involve coarse-graining of the charge distribution (51,52) or generating effective charges to reproduce the protein's external field (54). Given that we are interested in the electrostatic fluctuations, which are more sensitive to a local charge distribution, atomic charges from force-field potentials are more consistent with our purpose. Our approach starts with solving the standard ENM problem to obtain a set of displacements δr im for each normal mode m with the eigenvalue λ m . Each of these displacements is then applied to all charges q ik of residue i, where index k runs over all atoms of residue i (Fig. 1). Displacements δr im in turn produce dipole moments δµ (m) ik = q ik δr im which interact with charges of the protein active site. The overall electrostatic effect of the fluctuating protein medium is given as a sum over all contributions from each partial charge of all the residues except for the active site. This cumulative effect of all protein charges is given by the frequency-dependent response function of the electrostatic potential χ φ (ω) Here, E α 0i is the α Cartesian projection of the electrostatic field produced by all charges q ik of residue i at an atom in the active site ( Fig. 1) Here, r 0 is the position of the atom in the active site and r ik is the position of charge q ik of residue i.

Electrostatic field response function
The potential response function χ φ (ω) in Eq. (6) describes an oscillating electrostatic potential produced by the protein matrix in response to placing an oscillatory charge at the position in the active site where the potential is recorded (Fe ion in our MD simulations). Similarly, the response function of the electric field considers the field produced by the protein in response to an oscillating dipole moment m 0 (t) = m 0 (ω) exp[iωt] placed in the active site. The deformation of the protein matrix induced by this probe dipole results in the electric field E 0 (ω) at the same site. It can be found by summation over the contributions from all individual dipoles δµ ik (ω) arising from residues' atomic charges Here, T ik is the dipolar tensor connecting the position of the active site with the charge q ik of residue i. The dipole moment δµ ik (ω) at the position of q ik is caused by the residue displacement δr ik (ω) and can be expressed in terms of the displacement response function where δF j (ω) = ∑ k q jk T jk · m 0 (ω) is the force caused by the dipole m 0 (ω) at residue j. We therefore obtain a linear relation between the dipole m 0 (ω) placed at the active site and the electric field E 0 (ω) produced by the protein in response to this perturbation. The proportionality coefficient between the external perturbation and the response is the response function given by the rank-2 tensor Here, as above, the summation is done over the repeated Greek indices denoting the Cartesian components of the corresponding tensors. We will be mostly interested in the trace of the tensor

Response and correlation functions
The response functions calculated by the DENM formalism need to be related to time correlation functions supplied by the MD trajectories. The dynamics of a general dynamical variable δX (t) = X (t) − X is characterized by the time self-correlation function The normalized correlation function S X (t) = C X (t)/ (δX ) 2 is related to the response function by the fluctuation-dissipation equation which forms the basis of our analysis (55) whereS X (z) is the Laplace-Fourier transform of S(t) defined in the upper half of the complex plane of z.
Since our main focus is on variances of physical properties, we will define the generalized compliance for the variable X as follows The same property can be calculated from the ω = 0 value of the real part of the response function In case of residue displacements considered as variable X (X = r), λ r ∝ C −1 is proportional to the inverse force constant of the network springs. As for the compliance of a macroscopic body, defined as the inverse of stiffness, λ r depends on the shape and boundary conditions, in contrast to elastic moduli representing material properties. For the elastic network, the generalized compliance reports on how the softness of the network affects the variance of the variable of interest. In the case of X = φ reporting on the potential φ at the position of the heme's iron, λ = e 2 λ φ is the reorganization energy of a half redox reaction corresponding to changing the oxidation state of the heme (6,56). The standard definition of half-reaction reorganization energy involves, instead of one atom, the distribution of the electronic density of the transferred electron over a few atoms of the active site. We simplify this problem here by assuming all electron charge e localized at the centroid of this charge density, the iron atom of the heme. The electrostatic potential at a single atom is a well defined physical property, but its variance is only an approximate representation of the observable reorganization energy of changing the redox state (6).
We want to gain insight into the dissipative dynamics of the elastic network and its comparison with the dynamics calculated from MD simulations. The loss function χ ′′ X (ω) provides direct access to both the set of characteristic relaxation frequencies and their relative contributions. However, it does not weigh the low and high frequencies as they contribute to the variance in Eq. (14). Therefore, in addition to the loss function, we will consider the function This function is normalized, ∞ 0 α X (ω)dω = 1, and tends to the characteristic relaxation time τ in the ω → 0 limit

Proteins used in case studies and the simulation protocol
Three heme proteins, oxidized form of myoglobin (metmyoglobin, metMB, 1YMB) reduced state of cytochrome B562 (cytB, 256B), and reduced state of bovine heart cytochrome c (cytC, 2B4Z) were simulated by all-atom MD. The parameters of the proteins were taken from CHARMM27 force field (57), and NAMD (58) was used for the trajectories production. Each protein was placed at the center of a cubic box with the side length of ≃ 108Å and solvated with TIP3P waters (59). The number of waters used in simulations were: 32891 (metMB), 33268 (cytB), and 33189 (cytC). The VMD (60) plugin Autoionize was used to neutralize the simulation cell by adding Na + and Cl − ions to bring the ionic strength of the solution to 0.1. Particle-mesh Ewald with the grid resolution < 1Å was used for electrostatic interactions, and all other non-bonded interactions were calculated within 12Å cutoff. Following energy minimization, each protein was simulated for 5 ns in a NPT ensemble at P = 1 atm and T = 300 K. Temperature and pressure were controlled by using the Langevin dynamics with the damping coefficient of 5 ps −1 . The production NVE trajectory was started at the end of each 5 ns trajectory. The NVE ensemble is required for the proper sampling of the long-time dynamics since we found that using NPT and NVT ensembles artificially accelerates the observables relaxing on sub-nanosecond to nanosecond time-scale (61). The integration step was 2 fs and simulation frames for the analysis were saved each 0.05 ps. The simulation trajectories were 65 ns (metMB), 100 ns (cytC), and 123 ns (cytB). The force field standard parameters of the heme in the reduced state are taken from CHARMM27 (57). The charge of iron is q Fe = 0.24 in the reduced state. No standard parameters are available for the oxidized heme required for the simulations of metmyoglobin. These force field parameters were taken from Autenrieth et al (62). The atomic charge of oxidized Fe is q Fe = 1.34 in this parametrization. The iron atom is five-coordinate in its oxidized state. In order to allow the iron to shift out of the porphyrin plane, the harmonic force constant of the bond between Fe and N ε 2 of histidine 93 was adjusted to 65 kcal/mol as suggested in Ref. (63) and also implemented in Ref. (64).

Statistics and dynamics of residue displacements
Consistent with many previous studies, we have found that elastic networks are capable of reproducing the basic pattern of the distribution of mean-square displacements (msd's) along the protein backbone. The left panel in Fig. 2 compares residue msd's from MD with DENM calculations, while two other panels show the maps of variance-covariance matrices. Overall, there is a good agreement between the alteration in residue displacements along the backbone, and corresponding cross-correlations, calculated from the two sets of data. However, the network force constant required for the best fit of the msd from MD simulations, ≃ 0.3 kcal/(molÅ 2 ), is somewhat lower than the value of ≃ 1 kcal/(molÅ 2 ) typically adopted from fitting the Bfactors from crystallography (18)(19)(20). This value is also about a factor of two lower than 0.6 kcal/(molÅ 2 ) adopted below to globally fit the variances of the field and electrostatic potential fluctuations from MD data for all three proteins studied here.
The fitting of the force constant to the absolute msd's from MD makes the network too soft. This outcome is expected since the elastic network cuts high-frequency vibrations from the density of states and thus underestimates the absolute magnitudes of the residue displacements (27). It appears more consistent to use for the purpose of fitting only the portion of the msd related to protein's global motions. This component can in fact be separated from the vibrational part in the temperature dependence of the msd since global motions of the protein become observable at high temperatures, above the temperature of the protein dynamical transition (65). This high-temperature portion of the protein msd is about a half of the total magnitude (66), which is close to the factor required to reconcile the force constants from electrostatic variances and msd's. Figure 3 shows the dynamics of positions r i (t) of individual residues within the network and the dynamics of distances between residues r i j (t) = r i (t) − r j (t). It reports the loss functions χ ′′ i (ω) and χ ′′ i j (ω) calculated from the self-correlation functions of individual residues and self-correlation functions of distances between the residues These latter correlation functions are experimentally available from FRET (67) recording the evolution of the distance between two residues tagged with chromophores. The comparison of χ ′′ i (ω) from DENM and MD is illustrated in Fig. 3a for only one residue, i = 21, belonging to cytB protein. The pair of residues with i = 21 and j = 84 is used to calculate C i j (t) in Fig. 3b. The position of the pair in the backbone of cytB is shown in Fig. 4. Examples of calculations for several additional pairs of residues from cytB can be found in Supporting Materials.
Both correlation functions, C i (t) and C i j (t), report the dynamics on roughly two time-scales, seen as two peaks in the their loss functions in Fig. 3. The two peaks represent, with varying weights, the faster local dynamics of individual residues and the slower global dynamics of the network.
The elastic network with single-exponential, Debye dynamics does not capture two timescales of the dynamics from MD simulations. When the Debye memory kernelζ(ω) = ζ 0 is used in the response function χ m (ω) in Eqs. (4) and (5), the loss function shows only one relaxation peak, slightly deformed from a simple Debye form by the distribution of network's (i = 21) and C i j (t) (i, j = 21, 84, Eqs. (19) and (20)) and corresponding loss spectra χ ′′ i (ω) and χ ′′ i j (ω) shown in Fig. 3. The residues used in producing χ ′′ i j (ω) are connected by the red solid line. The red spheres representing the two chosen residues are centered at their corresponding C α carbons. The Fe atom of the heme is rendered as a yellow sphere. eigenvalues λ m . A single friction coefficient evidently misses the fact that energy dissipation decreases with increasing frequency of vibrations. Since the protein network is characterized by two spring constants, lower for non-covalent neighbors within the cutoff and higher for covalent neighbors, it must possess at least two characteristic friction coefficients.
Friction kernels typically used in applications are phenomenological (47) and we as well cannot offer a consistent theoretical formalism capturing the complex dynamics of residue displacements. We have therefore adopted a phenomenological response function best describing the MD results. Instead of searching for a functional form ofζ(ω) reproducing MD simulations, we have resorted to a two-Debye form for the entire response function of bead displacements in Eq. (4) Here, the amplitude a represents the relative weight of the high-frequency dissipation with the high-frequency friction ζ h . Correspondingly, ζ l is the low-frequency friction, which is larger in magnitude. It is obvious that Eq. (21) satisfies the static limit χ m (0) = λ −1 m . The two-Debye form of χ m (ω) is directly applied to calculate the response function of residue displacements where the function χ αβ i j (ω) is given by Eq. (5). The loss spectra of residue displacements show two peaks and qualitatively agree with the simulations. The positions and relative heights of the peaks can be adjusted by choosing ζ h,l and the amplitude a in Eq. (21), as is illustrated in Fig. 3a by dashed and dash-dotted lines. The model can therefore be well parameterized to reproduce the dynamics of displacements of a small subset of residues. However, it fails to discriminate between the dynamics of residues across the protein backbone.
The loss functions calculated from MD for a number of residues (see Supporting Materials) display similar two-peaks pattern, but the weights and positions of their peaks vary among the residues. In contrast, the dynamics of displacements of individual residues in the elastic network are mostly driven by the global motions of the entire network. As a result, there is little difference between the loss functions of individual residues calculated with DENM. The same statement applies to the dynamics of distances between residues shown in Fig. 3b. While one can reproduce the loss function of a chosen pair of residues by adjusting the parameters of χ m (ω) in Eq. (21), these parameters do not translate to all pairs in the network and instead can be only viewed as an average representation of the pairs dynamics.
The fast component of the dynamics, prominent in the loss function χ ′′ X (ω), is less significant in the function α X (ω) due to the 1/ω scaling of χ ′′ X (ω) in Eq. (16). The function α X (ω) is more relevant to problems related to the generalized compliance as defined by Eq. (14). The correct representation of the fast dynamics is therefore of lesser importance for these type of problems. The 1/ω scaling difference between χ ′′ X (ω) and α X (ω) explains the relative success of the network models in reproducing the pattern of residue msd's (Fig. 2). Those are obtained by frequency integration of α r (ω) and are less sensitive to the details of local dynamics of individual beads in the network represented by high-frequency peak of their loss functions.

Dynamics of the electrostatic field and potential
Essentially all time correlation functions obtained here from MD simulations show a pattern that is roughly represented by a two-exponential decay, requiring in some cases a third decaying exponent for a better mathematical fit. The relative weights of the fast and slow components of the correlation functions differ depending on the observable. The loss function of the electric field χ ′′ E (ω) at the position of protein's heme is mostly single-exponential, with the low-frequency peak much exceeding the high-frequency one. This pattern is less uniform for the loss function of the electrostatic potential showing different weights of the low-and high-frequency peaks in χ ′′ φ (ω) among the proteins (35). Figure 5 shows the loss spectra χ ′′ X (ω)/χ ′ X (0) calculated from the DENM and MD for cytB and metMB proteins. Similarly to the case with the residue displacements, the low-frequency coefficient needs adjustment to reproduce the position of the main loss peak: the value of ζ l = 30 ns was adopted for cytB and ζ l = 10 ns was taken for metMB. Still, the set of parameters taken to reproduce the loss spectrum of electrostatic potential does not perform as well when applied to the electric field (cf. Figs. 5c and 5d).
The relative contribution of the fast dynamics component is reduced in the spectral function α X (ω) describing the frequency variation of the generalized compliance (Eqs. (14 and (16)). The results for the electrostatic potential (α φ (ω), blue) and electric filed (α E (ω), red) are shown. A single set of network parameters is used for all three heme proteins: ζ l = 30 ns, ζ h = 0.006ζ l , a = 0.35, ε = 100. Figure 6 compares α X (ω) from DENM with MD data for all three globular proteins studied here. In order to show the ability of DENM to describe the entire set of MD data, a single set of network and friction parameters has been assigned to all three proteins. As expected, the performance of DENM is better for α X (ω). The average relaxation time τ , given by the ω → 0 limit of α X (ω) (Eq. (17)), is well reproduced by the network calculations. The deviations for some observables correspond to the cases when the dynamics are strongly dominated by its high-frequency component. For these cases, the fit can be improved by adjusting the weight a of the fast component in Eq. (21) to a higher value.
The loss functions χ X (ω) shown in Figs. 3, 5, and 6 are normalized with the corresponding values of χ ′ X (0) and therefore are not affected by the magnitude of the spring constant C describing the non-covalent interactions within the cutoff distance r c = 15Å. The absolute values of the variances are, however, proportional to k B T /C, as for instance in Eq. (2). We adopted in our calculations the value of k B T /C = 1Å 2 , which is equivalent to C ≃ 0.6 kcal/(molÅ 2 ). The Table 1: Generalized compliances (Eq. (14)) for the electrostatic potential fluctuations λ φ (also known as reorganization energies of redox half reaction λ = e 2 λ φ , kcal/mol) and for the electric force, F = eE (λ F , kcal/(molÅ 2 )) for the three proteins used in the case study; the network is assigned the force constant of C = 0.6 kcal/(molÅ 2 ) for non-covalent springs and εC, ε = 100 for covalently bound neighbors. The cutoff radius is 15Å. MetMB  85  228  436  243  CytB  115  136  9  17  CytC  141  159  18  83 results of calculating the generalized compliances λ X for electrostatic observables (Eqs. (14) and (15)) with this choice of the force constant are summarized in Table 1. The agreement is semi-quantitative at best. While the value of the force constant can clearly be adjusted in each particular case, this change does not propagate into equally good agreement between λ φ and λ F , where F = eE is the electrostatic force acting on the unit charge placed at Fe of the heme. The requirement to reproduce electrostatic forces might be more stringent than to describe the electrostatic potential because of a shorter range and higher sensitivity to the local structure of the former. Since the statistical variances listed in Table 1 are not affected by the modeling of dissipative dynamics, the properties of the network itself need to be further fine-tuned. The inclusion of softening of the residue displacements by electrostatic water fluctuations (14) appears to be the first avenue for the model improvement. This addition will allow one to accommodate for the effect of the interface, in addition to the motions dictated by the shape and mass distribution which elastic models capture in the first place (28)(29)(30). Incorporating the water response will effectively produce a more heterogeneous distribution of the network force constants.

CONCLUSIONS
Network models of biomolecules are coarse-grained representations designed to describe their large-scale collective conformational motions (27). They capture the basic topology and packing of residues in a biopolymer and robustly reproduce conformational global motions driven by the distribution of mass and shape (28)(29)(30). Electrostatic properties is yet another area where coarse-graining of biomolecules might be efficient. The long range of Coulomb interactions effectively averages out details of the local structure suggesting that one can potentially describe the statistics and long-time dynamics of electrostatic fluctuations by global motions of charges assigned to the elastic network. The demonstration that this approach is in principle consistent with all-atom MD simulations is the main result of this study. Two components are critical to our approach: force-field atomic charges distributed at the residues of the network and twoexponential, overdamped dynamics assigned to each normal mode diagonalizing the network Hessian.
Harmonic mechanical motions of the network of beads clearly do not incorporate dissipative dynamics of biomolecules in solution (34,45,46) and much is still needed to be done to achieve a physically consistent description of the protein dynamics. The elastic network employed here assigns weaker springs between all non-covalent neighbors within a cutoff radius and stronger springs between covalent neighbors. Correspondingly, two global friction coefficients are assigned to each normal mode diagonalizing the network Hessian. This phenomenological model qualitatively captures the two-peak loss spectrum of residue displacements and qualitatively similar loss spectra of the electrostatic potential and electric field. However, the characteristic adjustable friction coefficients used for the electrostatic fluctuations are not transferrable to the network displacements. A new set of parameters is needed when each property is considered separately.