Calculations

TB Thomas Barbot VB Veronica Beswick CM Cédric Montigny ÉQ Éric Quiniou NJ Nadège Jamin LM Liliane Mouawad

This protocol is extracted from research article:

Deciphering the Mechanism of Inhibition of SERCA1a by Sarcolipin Using Molecular Simulations

**
Front Mol Biosci**,
Feb 4, 2021;
DOI:
10.3389/fmolb.2020.606254

Deciphering the Mechanism of Inhibition of SERCA1a by Sarcolipin Using Molecular Simulations

Procedure

All the following calculations were carried out with program CHARMM (Brooks et al., 1983) using the CHARMM36 all-atom force field (Klauda et al., 2010; Best et al., 2012), excluding the CMAP correction, which is not needed for energy minimizations and normal mode calculations. The TIP3P model (Jorgensen et al., 1983; Neria et al., 1996) was used for water molecules. Electrostatic and van der Waals interactions were switched to zero with a cut-on distance of 6 Å and a cut-off distance of 10 Å, and the relative dielectric constant was equal to 3. These values were chosen to avoid important distortions of the structure due to extreme energy minimizations.

Each system (protein + membrane + water + ions, and protein + membrane in the absence of water) was energy minimized by 1,000 steps using the steepest descent method followed by conjugate gradient method for thousands of additional steps until the root mean square of the energy gradient (GRMS) fell below 10^{−4} kcal/mol.Å^{−1}.

The energy-minimized structures, E1.Mg^{2+}:SLN and E1.Mg^{2+}, were considered as starting systems to calculate normal modes. Because of the large number of atoms in each system (>70,000 atoms) the DIMB method (Diagonalization in a Mixed Basis) (Mouawad and Perahia, 1993; Perahia and Mouawad, 1995) was used in program CHARMM. DIMB is an iterative method in which the Hessian matrix to be diagonalized is divided into small-size Hessian matrices expressed in a mixed basis of intermediate NMs and Cartesian coordinates. By an iterative process of diagonalizations, the modes converge to those that would be obtained from the diagonalization of the entire Hessian matrix. However, to be able to perform the calculations in a reasonable time (<10 days), we accelerated this method by making it work on a GPU system, in a home-made adaptation of CHARMM. This adaptation was tested beforehand on small proteins to ensure that the results were the same as those obtained on CPU. For each system 206 modes were calculated, the 6 trivial global translation-rotation modes with null frequencies, and 200 non-trivial modes (modes 7 to 206), sorted according to the ascending order of their frequencies. The modes were considered as converged when the maximum eigenvectors convergence was <0.03.

The percentage of overlap (*p*) between two vectors, ${\overrightarrow{V}}_{1}$ and ${\overrightarrow{V}}_{2}$ was calculated using their dot product:

where $\Vert {\overrightarrow{V}}_{1}\Vert $ and $\Vert {\overrightarrow{V}}_{2}\Vert $ are the norms of vectors ${\overrightarrow{V}}_{1}$ and ${\overrightarrow{V}}_{2}$.

In this study, we calculated the percentage of overlap (or projection) between the modes of a structure that we call *A* and the difference coordinates vector between structures *A* and *B*. So ${\overrightarrow{V}}_{1}={\overrightarrow{q}}_{m}^{A}$, where ${\overrightarrow{q}}_{m}^{A}$ is mode *m* of structure *A*, with *m*∈[7, 206], and ${\overrightarrow{V}}_{2}=\Delta \overrightarrow{R}={\overrightarrow{R}}^{B}-{\overrightarrow{R}}^{A}$, where ${\overrightarrow{R}}^{A}$ and ${\overrightarrow{R}}^{B}$ are the coordinates of structures *A* and *B*, respectively. Vector $\Delta \overrightarrow{R}$ may also be written *A*→*B*. Prior to the projection, structure *B* is superimposed on structure *A* using C_{α} atoms. The projection is also done considering only C_{α} atoms in both vectors, ${\overrightarrow{q}}_{m}^{}$ and $\Delta \overrightarrow{R}$, in order to avoid misleading directions due to sidechains. The higher the percentage of overlap between these two vectors (${\overrightarrow{q}}_{m}$ and $\Delta \overrightarrow{R}$), the most straightforward is the transition from *A* to *B*. Here *A* refers to one of the two energy-minimized structures (E1.Mg^{2+} and E1.Mg^{2+}:SLN) and *B* to one of the two chosen crystal structures (1VFP for state E1.2Ca^{2+} and 3W5C for state E2). For the E1.2Ca^{2+} state, 1VFP was preferred to the other structures, 3AR2, 1T5S, and 3TLM, despite its slightly lower resolution in some cases, because 1T5S and 3TLM lack the disulfide bridge that stabilizes the protein, and the ATP analog in 3AR2 chelates Ca^{2+} instead of Mg^{2+} (see Supplementary Table 1). For the E2 state, the structure of 3W5C, although lacking an ATP analog, was preferred to 2DQS, because it is free of exogenous molecules like thapsigargin. 2DQS, which has the ATP analog, also contains thapsigargin to stabilize it and disrupt the communication between the transmembrane domain and the cytosolic headpiece of the ATPase (Picard et al., 2006; Montigny et al., 2007). However, despite these differences, the two E2 structures, 3W5C and 2DQS, are very similar, with an RMSD between all their C_{α} atoms of only 0.6 Å.

We also calculated the percentage of overlaps between all the modes calculated from one structure (*B*) over one chosen mode calculated from another structure (*A*), in which case, ${\overrightarrow{V}}_{1}={\overrightarrow{q}}_{m}^{A}$ and ${\overrightarrow{V}}_{2}={\overrightarrow{q}}_{n}^{B}$, where *m* = 54 when *A* = E1.Mg^{2+}:SLN and *n*∈[7, 206] for *B* = E1.Mg^{2+}, or *m* = 40 when *A* = E1.Mg^{2+} and *n*∈[7, 206] for *B* = E1.Mg^{2+}:SLN (see subsection The P-Domain Plays a Central Role in the Transition Toward the E2 State in the Results).

For each structure, E1.Mg^{2+}:SLN and E1.Mg^{2+}, the atomic fluctuations, *f*_{i}, were calculated from the 200 modes using:

where ${\overrightarrow{r}}_{i}$ is the position of atom *i*, ${\overrightarrow{q}}_{im}^{}$ its component in mode *m*, ω_{m} the frequency of mode *m, T* is the temperature, and *k*_{B} the Boltzmann constant. The sum runs from 7 to *l*, which is usually equal to 3*N*, where *N* is the number of atoms. However, the fluctuations are dominated by those of the lowest frequency modes, since ${f}_{i}^{2}$ is inversely proportional to the frequency. Here *l* = 206. But, even when all modes are considered, the fluctuations obtained from NM in the all-atom model are undervalued compared to those obtained from molecular dynamics or the crystal B-factors. Therefore, for this calculation, *T* was taken equal to 1,300 K to reach the range of the fluctuations obtained from the B-factors.

The fluctuations were also calculated from the crystal structure B-factors using the equation:

The correlation C_{ij} between two C_{α} atoms, *i* and *j*, calculated from one mode, *m*, is as follows:

All analyses were performed with CHARMM. The plots were drawn using Kaleidagraph (http://www.synergy.com) and the images were done using VMD (Visual Molecular Dynamics, https://www.ks.uiuc.edu/Research/vmd; Humphrey et al., 1996).

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Note: The content above has been extracted from a research article, so it may not display correctly.

Q&A

Your question will be posted on the Bio-101 website. We will send your questions to the authors of this protocol and Bio-protocol community members who are experienced with this method. you will be informed using the email address associated with your Bio-protocol account.