# Simulation studies of structural changes and relaxation processes in lysozyme under pressure

Véronique Hamon, Paolo Calligari, Konrad Hinsen, Gerald R. Kneller

► To cite this version:

Véronique Hamon, Paolo Calligari, Konrad Hinsen, Gerald R. Kneller. Simulation studies of structural changes and relaxation processes in lysozyme under pressure. Journal of Non-Crystalline Solids, 2006, 352 (42-49), pp.4417-4423. 10.1016/j.jnoncrysol.2006.01.141. hal-00388010

HAL Id: hal-00388010

https://hal.science/hal-00388010v1

Submitted on 25 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Simulation studies of structural changes and relaxation processes
in lysozyme under pressure

Véronique Hamon <sup>a</sup>, Paolo Calligari <sup>b,c</sup>, Konrad Hinsen <sup>b</sup>, Gerald R. Kneller <sup>a,b,*</sup>

<sup>a</sup> Centre de Biophysique Moléculaire, CNRS<sup>1</sup>, Rue Charles Sadron, 45071 Orléans, France

<sup>b</sup> Laboratoire Léon Brillouin, CEA Saclay, 91191 Gif-sur-Yvette, France

<sup>c</sup> Institut Laue-Langevin, 6 rue Jules Horowitz, BP 156, 38042 Grenoble, France

The paper describes preliminary results of a molecular dynamics simulation study on the influence of non-denaturing hydrostatic pressure on the structure and the relaxation dynamics of lysozyme. The overall compression and the structural changes are in agreement with results from recent nuclear magnetic resonance experiments. We find that moderate hydrostatic pressure reduces essentially the amplitudes of the atomic motions, but does not change the characteristics of the slow internal dynamics. The latter is well described by a fractional Ornstein-Uhlenbeck process, concerning both single particle and collective motions.

Keywords: Biopolymers; Neutron diffraction/scattering; Molecular dynamics; Structural relaxation

**1. Introduction**

Different experimental techniques, like quasielastic neutron scattering, fluorescence correlation spectroscopy (FCS), and kinetic studies of ligand rebinding show that protein dynamics is characterized by a vast spectrum of relaxation rates, ranging from picoseconds to hours [1–3]. The typical time scales observed by quasielastic neutron scattering are in the picosecond to nano-second range (10<sup>-12</sup>–10<sup>-9</sup> s), FCS probes protein dynamics in the millisecond to second range, and studying ligand rebinding after flash photolysis explores time scales in the millisecond to hour range. Xie and co-workers used the fractional Ornstein–Uhlenbeck process as a model to describe the relative motions of the flavin fluophore in flavin reductase with

respect to a fluorescence-quenching tyrosine side-chain in the vicinity [4]. Models describing so-called fractional Brownian dynamics are increasingly used to interpret the relaxation dynamics of complex systems which do not exhibit exponential relaxation, but a wide spectrum of relaxation rates [5].

Recently it has been shown that fractional Brownian Dynamics is also seen on the pico- to nano-second time scale, which is explored by quasielastic neutron scattering and molecular dynamics (MD) computer simulations [6,7]. Computer simulations play an important role in developing models for the dynamics of liquids and proteins, since they allow to perform numerical experiments on a complexity scale between analytical calculations and experiments. Using techniques known from signal processing, one can in particular compute memory functions of time correlation functions from the simulated trajectories [8]. The latter yield a rigorous description of relaxation processes on a microscopic basis [9].

Using the models and methods exposed in [8,6,7] we investigate here in a preliminary study if and how high,

<footnote n="1"><sup>1</sup> Corresponding author. Address: Centre de Biophysique Moléculaire, CNRS, Rue Charles Sadron, 45071 Orléans, France. Tel.: +33 2 38 25 78 42; fax: +33 2 38 63 15 17.</footnote>

<footnote>E-mail address: kneller@cnrs-orleans.fr (G.R. Kneller).</footnote>

<footnote n="1">¹ Affiliated with the University of Orléans.</footnote>
but non-denaturing hydrostatic pressure influences the internal structure and relaxation dynamics of lysozyme in solution. First, we present the theoretical background of relaxation phenomena in complex systems and then we show analyses of molecular dynamics (MD) simulations of a lysozyme protein solution at ambient temperature and different pressures. We compare the overall compression and the essential structural changes with results from recent nuclear magnetic resonance (NMR) experiments and investigate the influence of pressure on the relaxation dynamics. Here the single particle dynamics as well the collective dynamics are considered. The paper is concluded by a short discussion of the results.

## 2. Theoretical background

### 2.1. Memory function and relaxation time spectrum

In the following we consider an arbitrary dynamical variable $a$, whose time evolution can be described in terms of classical Hamiltonian dynamics. Such a dynamical variable depends on the momenta $p = \{p_1, ..., p_n\}$ and the coordinates $q = \{q_1, ..., q_n\}$ spanning the $2n$-dimensional phase space. Here $n$ is the number of degrees of freedom. Formally, the time dependence of $a$ is given by $a(t) = \exp(i\mathcal{L}t)a(0)$, where $\mathcal{L} = -i\left\{\frac{\partial H}{\partial p_i} - \frac{\partial H}{\partial q_i}\right\}$ is the Liouville operator and $H(p,q)$ is the Hamilton function of the system under consideration. Here $a(t)$ is a short notation for $a(p(t),q(t))$. We define the autocorrelation function

$$c(t) = \langle a^*(0)a(t) \rangle \quad (2.1)$$

as an ensemble average over the equilibrium distribution in phase space, where the asterisk denotes a complex conjugate. Zwanzig has shown that the time evolution of $c(t)$ is described by the integro-differential equation

$$\frac{dc(t)}{dt} = - \int_{0}^{t} d\tau \kappa(t - \tau) c(\tau), \quad (2.2)$$

where $\kappa(t)$ is the memory function which is associated to the time correlation function $c(t)$ [9,10]. Introducing the projector $\mathcal{P}$ on the variable $a$ of interest, whose action on another phase space variable $f \equiv f(p,q)$ is defined through

$$\mathcal{P}f = \frac{a\langle af \rangle}{\langle a^2 \rangle}, \quad (2.3)$$

$\kappa(t)$ can be expressed as

$$\kappa(t) = \frac{\langle \dot{a} \exp(i[1 - \mathcal{P}\mathcal{L}t]\dot{a}) \rangle}{\langle a^2 \rangle}. \quad (2.4)$$

The above expression shows that the memory function allows to describe relaxation processes, at least formally, on the basis of microscopic Hamiltonian dynamics.

For theoretical studies of slow relaxation processes it is useful to introduce the normalized autocorrelation function

$$\psi(t) = \frac{\langle a^*(0)a(t) \rangle}{\langle |a|^2 \rangle}, \quad (2.5)$$

writing the latter in the form

$$\psi(t) = \int_{0}^{\infty} d\lambda p(\lambda) \exp(-\lambda t). \quad (2.6)$$

Here $\lambda \ge 0$ are relaxation rates and $p(\lambda)$ is a normalized distribution with $\int_0^\infty d\lambda p(\lambda) = 1$. It follows from (2.6) that the Laplace transform of $\psi(t)$, which is defined as $\hat{\psi}(s) = \int_0^\infty dt \exp(-st)\psi(t)$ ($\Re\{s\} > 0$), is related to $p(\lambda)$ through a Stieltjes transform,

$$\hat{\psi}(s) = \int_{0}^{\infty} d\lambda \frac{p(\lambda)}{s + \lambda}. \quad (2.7)$$

The latter may be inverted to give [11]

$$p(\lambda) = \lim_{\epsilon \to 0^+} \frac{1}{\pi} \Im\{\hat{\psi}(-[\lambda + i\epsilon])\}. \quad (2.8)$$

It follows from the Laplace transform of the memory function equation (2.2) that $\hat{\psi}(s)$ has the simple form

$$\hat{\psi}(s) = \frac{1}{s + \hat{\kappa}(s)}, \quad (2.9)$$

and one obtains thus from (2.8) a relation between the spectrum of relaxation rates, $p(\lambda)$, and the memory function. It must be emphasized that the existence of a positive spectrum of relaxation times depends on the form of $\psi(t)$ and is thus not guaranteed. It is instructive to compute $p(\lambda)$ from (2.8) for the case of exponential relaxation, where $\psi(t) = \exp(-\lambda_0 t)$. As it should be, one finds $p(\lambda) = \delta(\lambda - \lambda_0)$, where $\delta(\cdot)$ denotes a Dirac distribution.

### 2.2. Fractional Ornstein–Uhlenbeck process

Any dynamical model for the internal dynamics of proteins must account for the fact that atomic motions in a protein are confined in space and that the relaxation dynamics in proteins is reflected in strongly non-exponential correlation functions. It has been shown by one of us (GRK) that the fractional Ornstein–Uhlenbeck process is an analytical model which describes the dynamics of single atoms in lysozyme, as well as its internal collective dynamics on the pico- to nano-second time scale [7] at normal pressure. The question we ask here is if and how the parameters of the model are influenced by applying an external pressure.

The classical Ornstein–Uhlenbeck process describes the diffusion of a particle in a harmonic potential ($K > 0$)

$$U(x) = \frac{1}{2} K x^2 \quad (2.10)$$

and leads to an exponentially decaying position autocorrelation function. In order to account for a broad spectrum of relaxation rates one considers a stochastic process which is described by a fractional Fokker–Planck equation,

$$\frac{\partial P(x, t)}{\partial t} = \tilde{\tau}^{1-\alpha} \mathcal{D}_t^{1-\alpha} \mathcal{L}_{FP} P(x, t). \quad (2.11)$$
Here $P(x,t) \equiv P(x,t|x_0,t_0)$ is the transition probability density for a move from $x_0$ at time $t_0$ to $x$ at time $t$, $\mathcal{L}_{\text{FP}}$ is the Fokker-Planck operator [12,13]

$$ \mathcal{L}_{\text{FP}} = \eta \frac{\partial}{\partial x} x + D \frac{\partial^2}{\partial x^2}, \qquad (2.12) $$

and $\mathcal{D}_t^{1-\alpha}$ denotes a fractional derivative of order $1-\alpha$ [14]. For an arbitrary function $f$ the latter is defined as

$$ \mathcal{D}_t^{1-\alpha} f(t) = \frac{d}{dt} \int_0^t d\tau \frac{(t-\tau)^{\alpha-1}}{\Gamma(\alpha)} f(\tau), \quad 0 < \alpha \le 1, \qquad (2.13) $$

where $\Gamma(\cdot)$ is the Gamma function [15]. One sees that $f$ is not altered in case of $\alpha = 1$. The two positive constants $\eta$ and $D$ in the Fokker-Planck operator (2.12) denote, respectively, a relaxation rate and the short-time diffusion constant. The parameter $\tilde{\tau}$ appearing in the right-hand side of (2.11) has the dimension of a time and must be introduced to assure its correct physical dimension.

The normalized position correlation function of the fractional Ornstein-Uhlenbeck process, $\psi(t) = \langle x(t)x(0) \rangle / \langle x^2 \rangle$, has the form

$$ \psi(t) = E_{\alpha}(-\eta_{\alpha} t^{\alpha}), \qquad (2.14) $$

where $E_{\alpha}(z)$ denotes the Mittag-Leffler function [16] and $\eta_{\alpha} = \tilde{\tau}^{1-\alpha}\eta$. The latter has the series representation

$$ E_{\alpha}(z) = \sum_{k=0}^{\infty} \frac{z^k}{\Gamma(1 + \alpha k)}. \qquad (2.15) $$

The exponential function is retrieved for $\alpha = 1$, where $\Gamma(1+k) = k!$. The memory function associated with $\psi(t)$ has the form [6]

$$ \kappa(t) = \frac{\alpha - 1}{\Gamma(\alpha)\tau_{\alpha}^2} \left(\frac{t}{\tau_{\alpha}}\right)^{\alpha-2}, \quad t > 0, \qquad (2.16) $$

for $0 < \alpha < 1$. Here $\tau_{\alpha} = \eta_{\alpha}^{-1/\alpha}$. Although $\kappa(t)$ is singular for $t \to 0+$, it is normalized to zero,

$$ \int_0^\infty dt \kappa(t) = 0. \qquad (2.17) $$

Such a behavior indicates that $\kappa(t)$ is a distribution – similar to the Dirac distribution – and not a normal function.

The Fourier spectrum of the correlation function and the associated spectrum of relaxation rates can be computed analytically [6,7]. One obtains a generalized Lorentzian for the Fourier spectrum,

$$ \tilde{\psi}(\omega) = \frac{2\tau_{\alpha} \sin(\alpha\pi/2)}{|\omega\tau_{\alpha}|(|\omega\tau_{\alpha}|^{\alpha} + 2\cos(\alpha\pi/2) + |\omega\tau_{\alpha}|^{-\alpha})}, \quad 0 < \alpha \le 1, \qquad (2.18) $$

and the following relaxation rate spectrum,

$$ p(\lambda) = \frac{\tau_{\alpha} (\tau_{\alpha} \lambda)^{\alpha-1} \sin(\pi \alpha)}{\pi (\tau_{\alpha} \lambda)^{2\alpha} + 2 (\tau_{\alpha} \lambda)^{\alpha} \cos(\pi \alpha) + 1}, \quad 0 < \alpha < 1. \qquad (2.19) $$

Both the time correlation function and its spectrum exhibit an algebraic behavior for large times and frequencies,

respectively. From the analytical properties of the Mittag-Leffler function it follows that [16]

$$ \psi(t) \propto t^{-\alpha} \qquad (2.20) $$

for $t \gg \tau_{\alpha}$, and expression (2.18) shows that

$$ \tilde{\psi}(\omega) \propto \omega^{-(1+\alpha)} \qquad (2.21) $$

for $\omega \gg \tau_{\alpha}^{-1}$. It should also be noted that

$$ p(\lambda) \propto \begin{cases} \lambda^{\alpha-1} & \text{if } \lambda \ll \tau_{\alpha}^{-1}, \\ \lambda^{-(1+\alpha)} & \text{if } \lambda \gg \tau_{\alpha}^{-1}. \end{cases} \qquad (2.22) $$

The first of the above relations shows that the relaxation spectrum is not bound for relaxation rates tending to zero, which correspond to relaxation times tending to infinity.

Fig. 1 shows the function $\psi(t)$ given by Eq. (2.14) for $\alpha = 0.5$ (solid line), as compared to a normal exponential

<figure><img src="image_2.png" /><figcaption>Fig. 1. Model correlation functions $\psi(t) = E_{\alpha}(-t^\alpha)$ for $\alpha = 0.5$ (solid line), $\psi(t) = \exp(-t)$ (broken line), and $\psi(t) = \exp(-t^\alpha)$ for $\alpha = 0.5$ (dashed-dotted line). The inset shows the memory function corresponding to $\psi(t) = E_{\alpha}(-t^\alpha)$ for $\alpha = 0.5$. More explanations are given in the text.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Relaxation rate spectrum (2.19) for $\alpha = 0.5$ (full line) and $\alpha = 0.98$ (broken line).</figcaption></figure>
function (broken line) and a stretched exponential function, $\exp(-t^\alpha)$, for $\alpha = 0.5$ (dashed-dotted line). The insert displays the memory function (2.16) which is associated with the correlation function (2.14). Here the memory function is given by expression (2.16) for $t > \epsilon = 0.1$, and for $0 \le t < \epsilon$ it is represented by a linear function, in order to fulfil the normalisation condition (2.17). In the limit $\epsilon \to 0$, where $\kappa(t)$ becomes a distribution, the concrete representation in the interval $0 \le t < \epsilon$ does not matter. Fig. 2 shows the relaxation rate spectrum corresponding to the correlation function (2.14) for $\alpha = 0.5$ and for $\alpha = 0.98$. In the latter case one recognizes that $p(\lambda)$ is centered on $\lambda = \tau_\alpha^{-1}$.

### 3. Simulation and analysis

#### 3.1. Simulation protocol

We have simulated a single lysozyme molecule in water, using an orthorhombic simulation box of dimensions $6.16 \times 4.19 \times 4.61$ nm³. The coordinates of lysozyme have been taken from entry 193L [17] of the Brookhaven Protein Data Bank [18]. From this coordinate set the positions of the hydrogen atoms were reconstructed, leading to 1960 atoms for the lysozyme molecule. The simulation box was filled with 3403 water molecules, containing thus 12169 atoms in total.

All MD simulations were performed with the MMTK simulation package [19], using the Amber94 force field [20]. In order to account for realistic thermodynamic boundary conditions, all simulations were performed in the thermodynamic NpT-ensemble, with fixed number of particles and prescribed averages for temperature and pressure [21,22]. In all cases the temperature was chosen to be $T = 300$ K and the pressures were varied from ambient pressure (1 bar) up to 7 kbar. Here only simulations at 1 bar, 2 kbar, and 3 kbar are reported on. The simulated trajectories were recorded for 1.2 ns at $p = 1$ bar and for 1 ns at $p = 2$ kbar and $p = 3$ kbar, using in all cases a sampling time step of $\Delta t = 0.04$ ps for later analyses. All simulations have been carried out with an integration time step of 1 fs.

#### 3.2. Trajectory analysis

In order to compare structural changes with results obtained by structural NMR, which have been reported in the literature [23], we computed a difference distance matrix from the average structure of lysozyme at ambient pressure and at $p = 2$ kbar. These calculations have been performed using basic modules from the MMTK package.

The analysis of the relaxation dynamics on the basis of the recorded MD trajectories has been performed with the package nMoldyn [24]. Here we compare simulations at ambient pressure and at $p = 3$ kbar. In order to investigate the influence of pressure on the single particle relaxation dynamics we computed the average atomic mean-square displacement

$$W(t) = \frac{1}{N} \sum_{j=1}^{N} w_j \langle [\mathbf{R}_j(t) - \mathbf{R}_j(0)]^2 \rangle \quad (3.23)$$

at different pressures. Here $N$ is the number of atoms under consideration, and $w_j$ are positive weights, with $\sum_{j=1}^N w_j = N$. We chose $w_j$ to be proportional to the squared incoherent neutron scattering lengths [25], which means in practice to consider essentially the hydrogen atoms. We concentrate on the hydrogen atoms since their motions have on average the largest amplitudes and since the motion of hydrogen atoms is probed in quasielastic neutron scattering experiments, which will be used later for comparison in a more detailed study.

Corresponding changes in the collective relaxation dynamics have been examined by studying the coherent intermediate scattering function

$$c(\mathbf{q}, t) = \langle \delta\rho^*(\mathbf{q}, 0)\delta\rho(\mathbf{q}, t) \rangle, \quad (3.24)$$

which is computed from the spatially Fourier transformed atomic density fluctuation

$$\delta\rho(\mathbf{q}, t) = \rho(\mathbf{q}, t) - \langle \rho(\mathbf{q}, t) \rangle. \quad (3.25)$$

Here $\rho(\mathbf{q}, t) = \sum_{j=1}^N w_j \exp(i\mathbf{q} \cdot \mathbf{R}_j(t))$ is the Fourier-transformed atomic density and $\mathbf{R}_j(t)$ are the positions of the atoms in the lysozyme protein as a function of time. Using $\delta\rho(\mathbf{q}, t)$ and not $\rho(\mathbf{q}, t)$ as dynamical variable avoids a large elastic peak in the Fourier transform of $c(\mathbf{q}, t)$ with respect to time, which is due to static correlations in the atomic positions. The weights $w_j$ have been chosen proportional to the coherent scattering lengths for neutron scattering [25]. The solvent molecules are not considered in the analysis. The correlation function $c(\mathbf{q}, t)$, its Fourier transform, and the associated memory function are analysed in terms of an autoregressive (AR) model for the dynamical variable $\delta\rho(\mathbf{q}, t)$,

$$\delta\rho(\mathbf{q}, t) = \sum_{k=1}^{P} a_k(\mathbf{q}) \delta\rho(\mathbf{q}, t - k\Delta t) + \epsilon(\mathbf{q}, t). \quad (3.26)$$

Here $P$ is the order of the AR process, $a_k(\mathbf{q})$ are $\mathbf{q}$-dependent coefficients, and $\epsilon(\mathbf{q}, t)$ is white noise with zero mean and variance $\sigma^2(\mathbf{q})$. The set of $P+1$ coefficients $\{a_k, \sigma^2\}$ is obtained from a fit to the MD trajectory of $\delta\rho(\mathbf{q}, t)$. Details are described in [8], and the method is implemented in the nMoldyn package [24]. A summary can be found in [7]. Here we note only two key features of the method:

1. AR time series modelling allows us to compute reliable estimates for memory functions by solving the discrete memory function equation

$$\frac{c(\mathbf{q}, n+1) - c(\mathbf{q}, n)}{\Delta t} = - \sum_{k=0}^{n} \Delta t \kappa(\mathbf{q}, n-k) c(\mathbf{q}, k) \quad (3.27)$$

for $\kappa(\mathbf{q}, n) \equiv \kappa(\mathbf{q}, n\Delta t)$.

2. AR time series modelling allows us to obtain high resolution spectra directly from the parameters $\{a_k, \sigma^2\}$, without applying numerical Fourier transforms and
window functions to reduce the numerical error of the estimated spectra. Defining $S(\mathbf{q}, \omega)$ to be the Fourier transform of $c(\mathbf{q}, t)$, we have

$$S(\mathbf{q}, \omega) = \frac{\Delta t \sigma^2(\mathbf{q})}{(1 - \sum_{k=1}^{P} a_k(\mathbf{q}) \exp[-i\omega k \Delta t]) (1 - \sum_{m=1}^{P} a_m^*(\mathbf{q}) \exp[i\omega m \Delta t])} \quad (3.28)$$

In the literature this approach is well known as the maximum entropy estimation of power spectra [26,27].

In the analyses shown in the following we used an AR model of order $P = 1000$.

# 4. Results

## 4.1. Compression and structural changes

We will first briefly report on the overall compression and on structural changes in lysozyme under pressure. Lysozyme is a protein whose structure has been intensively studied in the past. Of particular importance for our work are studies of lysozyme under pressure by X-ray crystallography [28] and by NMR [29,23]. The NMR studies are particularly useful for comparison with MD simulations since they are performed on proteins in solution. From our simulation studies we find for the simulated lysozyme molecule at $p = 3$ kbar a volume reduction of $\Delta V \approx -2\%$ with respect to ambient pressure [30]. Refaee et al. give an compressibility of 7.5/Mbar for the whole protein [23], leading thus to $\Delta V \approx -2.25\%$ at $p = 3$ kbar. Our simulation result is thus in excellent agreement with the NMR experiments. In [29] Katamari et al. state a shortening of the hydrogen bonds between the amid groups of the backbone of 0.02 Å on average at $p = 2$ kbar. This value is close to the average reduction of 0.03 Å we find from our simulations.

More detailed information on the influence of pressure on the structure of lysozyme can be obtained from the difference distance map for the $C_\alpha$-atoms in lysozyme which is given by Refaee et al. [23]. Fig. 3 shows the difference distance plot displayed in Ref. [23] (left part) together with our results from MD simulation (right part). In both cases the structure at $p = 2$ kbar is compared to the one at ambient pressure. In contrast to the NMR experiment, where two average structures are compared, which have been averaged over a small ensemble (50 structures), we compare average structures for whole MD trajectories. We note here that the NMR structures are deposited as entries 1GXX (2 kbar) and 1GXV (ambient pressure) in the Brookhaven Protein Databank. The coloring scheme is chosen such that red<fnref n="2" /> indicates a reduction of distances under pressure and blue an increase. The coloring scheme

in the secondary structure of lysozyme (upper left part of the figure) is the same as in the difference distance plot. Note that the simulation results are displayed in form of the full symmetric difference distance matrix. It is difficult to compare simulation and experiment on an absolute scale, since the structures to be compared have not been obtained in the same way, but the qualitative agreement is quite satisfactory. The comparison shows in particular in both cases residue 46, which is in a so-called β-turn in the more flexible ‘β-domain’ of lysozyme, which contains essentially loops and β-strands as secondary structure motifs, is moved closer to the rest of the structure. A similar observation obtained from the NMR results for residue 85, which is located in a short β-strand between two α-helices is less well visible in the MD results.

## 4.2. Single particle dynamics

Fig. 4 shows the average atomic mean square displacements in lysozyme at pressures of $p = 1$ bar (diamonds) and $p = 3$ kbar (circles). Assuming the model of fractional Brownian diffusion in a harmonic well, we obtain the analytical expression

$$W(t) = 2\langle\mathbf{u}^2\rangle(1 - E_\alpha(-[t/\tau_\alpha]^2)) \quad (4.29)$$

for the mean square displacement. Here it has been used that for motions confined in space $W(t) = 2\langle\mathbf{u}^2\rangle(1 - \psi(t))$. The horizontal lines in Fig. 4 indicate the plateau values of $W(\infty) = 0.037 \text{ nm}^2$ for $p = 1$ bar (solid line) and $W(\infty) = 0.026 \text{ nm}^2$ for $p = 3$ kbar (broken line). These values have been obtained by estimating the corresponding average atomic fluctuations by independent calculations from the MD trajectories. The model function (4.29) has thus two free parameters, $\alpha$ and $\tau_\alpha$, which were fitted in the statistically safe region $0 < t < 50$ ps. The resulting fits gave $\alpha = 0.51$ and $\tau_\alpha = 32.03$ ps at $p = 1$ bar (solid line) and $\alpha = 0.55$ and $\tau_\alpha = 29.41$ ps at $p = 3$ kbar (broken line). The results show that pressures up to 3 kbar lead essentially to a reduction of the amplitudes of the atomic motions in lysozyme, but not to an essential change in the internal dynamics.

A remark concerning the fit of the model (4.29) to the data obtained from MD simulation is in place here. It must be emphasized that $W(t)$ is a dynamical quantity, in contrast to the atomic fluctuation $\langle\mathbf{u}^2\rangle$. Suppose the MD trajectory is long enough to sample the atomic fluctuations reasonably well, such that the error is below 10%. The mean square displacement attains 90% of its plateau value $W(\infty)$ after $t \approx 1$ ns and 95% after $t \approx 4$ ns if we use the model (4.29) for $\alpha = 1/2$. To compute $W(t)$ with sufficient statistical accuracy for such times, the corresponding MD trajectories must be at least an order of magnitude longer. Given that most MD trajectories for protein simulations do not exceed a few ns, these considerations show that a precise calculation of mean square displacements of atoms in proteins for times beyond 100–200 ps is usually not feasible. The quasi-linear evolution of the mean square

<footnote n="2">² For interpretation of color in Fig. 3, the reader is referred to the web version of this article.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Left: Secondary structure of lysozyme and difference distance plot for the C<sub>α</sub>-atoms at ambient pressure and at *p* = 2 kbar given in Ref. [23] (with permission from Elsevier). The coloring scheme uses red for a reduction of the distance under pressure and blue for an increase. Right: The corresponding difference distance plot obtained from the MD structures averaged over the respective trajectories.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Mean square displacement of lysozyme in solution obtained from MD simulation at *p* = 1 bar (diamonds) and at *p* = 3 kbar (circles). The fitted solid and broken line display, respectively, the mean square displacement according to the model of fractional Brownian dynamics for *p* = 1 bar (*α* = 0.51, τ<sub>α</sub> = 32.03 ps) and *p* = 3 kbar (*α* = 0.55, τ<sub>α</sub> = 29.41 ps). The corresponding horizontal lines indicate the plateau values which have been estimated by a separate calculation from the MD trajectories.</figcaption></figure>

displacements depicted in Fig. 4 for *t* > 100 ps illustrates
this point. Such a behavior indicates free diffusion and is
clearly an artefact, since the motions of atoms in a protein
are confined.

4.3. *Collective dynamics*

In Fig. 5 we show the results for the coherent dynamic
structure factor for lysozyme in solution at *p* = 1 bar (thin

<figure><img src="image_6.png" /><figcaption>Fig. 5. Fourier spectrum of the autocorrelation function *c*(*q*, *t*) for the atomic density fluctuation of lysozyme at |*q*| = 10 nm<sup>-1</sup>. The thin black line corresponds to *p* = 1 bar and the thick grey line to *p* = 3 kbar. The broken line corresponds to a fit of the model spectrum (2.18) corresponding to the model of fractional Brownian motion for *δρ*(*q*, *t*), and the inset shows the corresponding memory functions.</figcaption></figure>

black line) and *p* = 3 kbar (thick grey line). The coherent
dynamic structure factor is here defined as Fourier trans-
form of *c*(*q*, *t*) with respect to time. To interpret the data
in terms of an analytical model we use again the fractional
Ornstein–Uhlenbeck process, except that the dynamical
variable is here the Fourier transformed atomic density,
i.e. a dynamical variable describing collective motions,
and not the position of a single atom. Using this model,
we assume that there is an effective force proportional to
*δρ* which drives the latter back to zero, i.e. to a homoge-
nous atomic density in space. Only the data corresponding
to *p* = 1 bar have been used for the model fit of the
dynamic structure factor and the corresponding memory
function (broken lines), and the corresponding data have
been taken from Ref. [6]. Fig. 5 shows that the trends of
the dynamic structure factors at *p* = 1 bar and *p* = 3 kbar
and the corresponding memory functions are very similar.
As for the case of single particle dynamics, we do not find a
noticeable influence of the exerted pressure on the collec-
tive relaxation dynamics of lysozyme. At this point it
should be mentioned that the AR model gives reliable esti-
mates for correlation functions and the associated memory
functions, which would allow to detect small changes in the
characteristics of the dynamics of the system under consid-
eration. Although it is difficult to associate explicit error
bars, one can perform analyses with different numbers of
poles and different sampling steps, which show consistent
results with differences in the percent range [8,7].

5. Conclusion

The preliminary analysis of molecular dynamics simula-
tions of lysozyme in solution under hydrostatic pressure
showed that moderate, non-denaturing pressure leads to
an overall compression which is in excellent agreement with
results obtained from structural NMR experiments. Differ-
ence distance plots for the Cα-atoms in the protein main
chain at *p* = 2 kbar show satisfactory agreement with cor-
responding results from NMR. These findings are encour-
aging to study the influence of pressure on the relaxation
dynamics of lysozyme by MD simulation. Analysing simu-
lations at ambient pressure and at *p* = 3 kbar we find that
moderate pressures below the denaturation threshold
reduce the atomic fluctuations by about 30%, but do not
influence the relaxation dynamics of the protein. This con-
cerns the single particle dynamics, which is reflected in the
time evolution of the mean square displacement, as well as
the collective dynamics. Both types of dynamics are well
described by a fractional Ornstein-Uhlenbeck process,
where the atomic position is the dynamical variable in
the case of particle dynamics and the atomic density fluctu-
ation in the case of collective dynamics. Our findings are
consistent with a preliminary analysis of data obtained
from quasielastic neutron scattering on lysozyme in solu-
tion, which were obtained from the IN5 spectrometer at
the Institut Laue-Langevin in Grenoble [31].

References

[1] W. Doster, S. Cusack, W. Petry, Nature 337 (1989) 754.

[2] H. Yang, G. Luo, P. Karnchanaphanurach, T. Louie, I. Rech, S. Cova, L. Xun, X. Xie, Science 302 (2003) 262.

[3] W. Glöckle, T. Nonnenmacher, Biophys. J. 68 (1995) 46.

[4] H. Yang, X. Xie, J. Chem. Phys. 117 (2002) 10965.

[5] R. Metzler, J. Klafter, Phys. Rep. 339 (2000) 1.

[6] G. Kneller, K. Hinsen, J. Chem. Phys. 121 (2004) 10278.

[7] G. Kneller, Phys. Chem. Chem. Phys. 7 (2005) 2641.

[8] G. Kneller, K. Hinsen, J. Chem. Phys. 115 (2001) 11097.

[9] R. Zwanzig, Statistical mechanics of irreversibility, Lectures in Theoretical Physics, Wiley-Interscience, New York, 1961.

[10] J. Boon, S. Yip, Molecular Hydrodynamics, McGraw Hill, 1980, see eqs. (3.5.23), (3.5.33), (3.5.37) and the corresponding references.

[11] J. Schwarz, J. Math. Phys. 46 (2005) 013501.

[12] C. Gardiner, Handbook of Stochastic Methods, 2nd Ed., Springer Series in Synergetics, Springer, Berlin, Heidelberg, New York, 1985.

[13] H. Risken, The Fokker-Planck Equation, second ed., Springer Series in Synergetics, Springer, Berlin, Heidelberg, New York, 1996.

[14] K. Oldham, J. Spanier, The Fractional Calculus, Academic, New York, London, 1974.

[15] M. Abramowitz, I. Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1972.

[16] A. Erdélyi, W. Magnus, F. Oberhettinger, F. Tricomi, Higher Transcendental Functions, McGraw Hill, 1955.

[17] M. Vaney, S. Maignan, M. RiesKautt, A. Ducruix, Acta Cryst. D Biol. Cryst. 52 (1996) 505.

[18] H. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. Bhat, H. Weissig, I. Shindyalov, P. Bourne, Nucleic Acids Res. 28 (2000) 235.

[19] K. Hinsen, G. Kneller, Mol. Sim. 23 (2000) 275.

[20] W. Cornell, P. Cieplak, C. Bayly, I.I.R. Gould, K. Merz Jr., D. Ferguson, D. Spellmeyer, T. Fox, J. Caldwell, P. Kollman, J. Am. Chem. Soc 117 (1995) 5179.

[21] H. Andersen, J. Chem. Phys. 72 (1980) 2384.

[22] S. Nosé, J. Chem. Phys. 81 (1984) 511.

[23] M. Refaee, T. Tezuka, K. Akasaka, M. Williamson, J. Mol. Biol. 327 (2003) 857.

[24] T. Rog, K. Murzyn, K. Hinsen, G. Kneller, J. Comp. Chem. 24 (2003) 657.

[25] S. Lovesey, Theory of Neutron Scattering from Condensed Matter, vol. I, Clarendon, Oxford, 1984.

[26] J. Burg, PhD thesis, Stanford University, Stanford (CA), USA, 1975.

[27] A. Papoulis, Probablity Random Variables Stochastic Processes, third ed., McGraw Hill, 1991.

[28] C. Kundrot, F. Richards, J. Mol. Biol. 193 (1987) 157.

[29] Y. Kamatari, H. Yamada, K. Akasaka, J. Jones, C. Dobson, L. Smith, Eur. J. Biochem. 268 (2001) 1782.

[30] V. Hamon, PhD thesis, Université d'Orléans, 2004.

[31] G. Kneller, M.-C. Bellissent-Funel, V. Hamon, Tech. Rep. 8-04-286, Institut Laue-Langevin, 2004.