WAXSiS: a web server for the calculation of SAXS/WAXS curves based on explicit-solvent molecular dynamics

SWAXS calculation

The WAXSiS server computes the SWAXS curves following the methods described in a recent article (9). Accordingly, WAXSiS runs an explicit-solvent MD simulation of the biomolecule, typically for 20–500 ps depending on the size of the biomolecule. During the simulation, position-restraining potentials (force constant: 1000 kJ/mol nm²) are applied to the backbone atoms of the biomolecule and on heavy atoms of ligands. That procedure ensures that the simulation samples conformations close to the initial structure, while allowing thermal fluctuations of side chains, water and counter ions. The scattering contribution from the excluded solvent is computed from an MD trajectory of a pure-water simulation system, which is stored on the WAXSiS server.

After the solute simulation has finished, a spatial envelope is constructed that encloses the solute at a preselected distance (7 Å by default). Our algorithm constructs the envelope from an icosphere, which is obtained from a regular icosahedron by recursively subdividing its triangluar faces into four smaller triangles. An icosphere after four such recursions is shown in Figure 1A. The envelope is then obtained by moving the vertices of the icosphere in a radial direction, until all vertices have a distance of at least d from all solute atoms in all simulation frames. Once constructed, the same envelope is applied throughout the remaining calculation. Molecular images in Figure 1 were rendered with PyMol (19).

The electron density of each simulation frame is decomposed into density inside and outside of the envelope:

(3)

(4)

and denote the electron density in the solute and pure-water simulation systems, respectively, and the subscripts i and o indicate density inside and outside of the envelope. An example for the atoms contributing to and is shown in Figure 1B. As shown previously (9,17), the net intensity is then given by

(5)

(6)

Here, and are the Fourier transforms of the densities and , respectively. 〈·?^(ω) denotes the average over the simulation frames, while the superscript highlights that this average is taken at a fixed solute orientation ω. The asterisk indicates complex conjugate. denotes the orientational average, taking account of the fact that the solutes are randomly oriented in solution. Hence, Equations (5) and (6) show that only the densities inside of the envelope contribute, justifying why I(q) can be computed from a finite simulation system. The three terms in Equation (6) have an intuitive interpretation; the first term corresponds to the scattering of the solute and its hydration layer, and the second to the scattering of the excluded solvent. The third term corresponds to correlations between bulk water and the density contrast inside the envelope (9). If the buffer subtraction scheme of Equation (2) is applied, I(q) is corrected following the procedure in (18).

Given the coordinates of atoms within the envelope, the scattering amplitude for an individual simulation frame is written as

(7)

where N_A is the number of atoms within the envelope, f_j(q) are the atomic form factors, and is the coordinate of atom j. The f_j(q) are computed following , where a_k, b_k, c are the Cromer–Mann parameters (20), which are published in tables (21). The analogous relation is applied to compute using the atoms inside the envelope of the pure-water system (see Figure 1B). In order to account for electron-withdrawing effects in water molecules, we apply the correction proposed by Sorenson et al. to the form factors of water atoms (22).

The orientational average is evaluated numerically. For each absolute value of the scattering angle q, a set of vectors (j = 1, …, J_q) is distributed uniformly on the surface of a sphere with radius q. The vectors are distributed following the spiral method, as done previously (9,17). In our previous article we tested the convergence of I(q) with J_q, using various q and different solutes. In line with theoretical considerations (23), we observed that J_q should be taken proportional to (qD)² to achieve a preselected convergence, where D denotes the maximum diameter of the solute. In WAXSiS, we therefore use J_q = Max{100, 0.2 · (qD)²}, leading to small uncertainties over the entire q-range.

The computed SWAXS curves are at small angles very sensitive to small alternations in the solvent density. Because the density of the experimental buffer might deviate from the density of the explicit solvent model, we correct the solvent density following the method described previously (9). In short, we add a small uniform electron density to the water density, both in the bulk water and in the solute simulation system, to match the bulk density with a preselected buffer density. By default, we use a buffer density of 334 e nm⁻³, corresponding to pure water, but the WAXSiS user may specify a different buffer density. It is important to note that we do not fit the solvent density to match the computed with the experimental intensity.

Fitting procedure

If the user uploads an experimental SWAXS curve I_exp, it is fitted to the predicted curve (and not vice versa), because the predicted curve does not contain any free parameters. WAXSiS fits I_exp by minimizing two different metrics,

(8)

(9)

where N denotes the number of q-points, σ_exp(q_i) are the experimental errors, and f and c the fitted parameters. Thus, besides the overall scale f, only a constant offset c is fitted to absorb some uncertainty in the buffer subtraction. The first metric (Equation 8) is weighted by the experimental errors and is widely used in SWAXS analysis, but it may impose spuriously high weights to small angles. The second metric (Equation 9) imposes more uniform weights over small and wide angles, and we found it to be particularly suitbable to interpret wide-angle data. The fitted curve and χ-value for each metric are reported by WAXSiS.

MD simulation setup

The MD simulations are conducted by the YASARA Dynamics software (YASARA Biosciences, Vienna, Austria). YASARA is capable of setting up MD simulations automatically from PDB structures, thereby adding missing hydrogen and heavy atoms, fixing steric problems, and correcting for other frequent imperfections in PDB files. Protein and nucleic acids are described by the AMBER03 force field (24), and water parameters are taken from the TIP3P model (25). Force field parameters for modified side chains or ligands are derived by YASARA's AutoSMILES method. Should the structure contain exotic metal ions for which no force field parameters are available, WAXSiS replaces them with Fe^{2 +} before running YASARA. For the calculation of SWAXS curves, however, the form factors of the original exotic ions are applied. Electrostatic interactions are computed using the particle-mesh Ewald method, and the dispersive interactions are described by a Lennard-Jones potential with a cutoff at 9 Å. Temperature is controlled at 298.15 K and the pressure at 1 bar using YASARA's default temperature and pressure control settings. Simulation frames are written every 0.5 ps to ensure that the solvent configurations are reasonably uncorrelated. Based on the convergence assessment in previous work (9), we compute the total number of generated MD frames as . Note that N_fr can also be controlled by the user option “Convergence” (see below). Before collecting simulation frames for the SWAXS calculation, the simulation system is equilibrated for 5% of the total simulation time (but not shorter than 3 ps).

Specifying the initial structure

The user may either upload structure of a biomolecule or specify a protein data bank (PDB) code. In the latter case, WAXSiS picks the biological assembly of that PDB code from the PDB. If multiple biological assemblies are stored for this PDB entry, the assembly may be specified. If no biological assembly is present, the user is notified and the crystallographic structure is taken instead. Finally, if the PDB code corresponds to an NMR ensemble, the first structure from the ensemble is taken, leading to a note in the job results.

Many PDB structures contain crystallization agents, or methionine residues might have been replaced with selenomethionine to aid solving the structure. Because such modifications are typically not present in the SWAXS sample, WAXSiS by default replaces selenomethionine with methionine side chains. In addition, WAXSiS tries to identify crystallization agents by default and removes them. Because PDB files do not contain reliable information on whether a molecule is a crystallization agent or a biologically relevant ligand, WAXSiS tries to detect crystallization agents using a combination of ligand size, number of atomic contacts, and residue name. Specifically, molecules are considered biologically relevant if one of the following criteria matches: (i) they form a covalent bond the biomolecule; (ii) they form more than 2.5 contacts per heavy atom to the biomolecule (contact cutoff 3.5 Å); (iii) they contain >20 heavy atoms; (iv) the residue name is HEM, indicating a heme ligand. Molecules with residue names GOL, BNG, BOG, PG, PE, PGE typically indicate glycerol, B-nonylglucoside, B-octylglucoside, and polyethylene glycol, which are always considered as crystallization agents. These criteria worked reliably for many test cases. Both the handling of crystallization agents and selenomethionine can be specified by the user. If in doubt, we recommend that the user removes unwanted molecules manually, uploads the structure, and specifies WAXSiS to 'Keep both ligands and crystallization agents' in the user options.

Input options

The following basic options can be specified by the user. Default values are printed in italic font:

• Email address (optional): if an email address is entered, the user will be notified when the job has finished.

• Ligand handling: (i) keep ligands, try to remove crystallization agents, (ii) keep ligands and crystallization agents, (iii) remove all.

• Buffer subtraction method: (i) buffer scattering reduced by solute volume (Equation 2), (ii) total buffer scattering subtracted (Equation 1).

• Maximum scattering vector q, defined as q = 4πλ⁻¹sin θ, where 2θ is the scattering angle (default: 1 Å⁻¹)

• Optional upload of an experimental SWAXS curve. The q units (Å⁻¹ or nm⁻¹) and the scattering convention (q or s = 2λ⁻¹sin θ) adopted for the experimental data can be specified.

The following advanced options can be specified:

• Output q units (Å⁻¹ or nm)

• Electron density of the buffer (default: 334 e nm⁻³)

• Replace selenomethionine with methionine (yes/no)

• Distance of envelope from the solute (default: 7 Å)

• Accuracy of convergence (quick/normal/thorough). Selecting 'quick' or 'thorough' causes WAXSiS to run the MD simulation four times shorter or five times longer as compared to 'normal' convergence, respectively. “Thorough” jobs run at lower priority in the queue.

• Generate a new random seed for initial velocities. Choosing yes/no leads to slightly varying or reproducible results, respectively.

Output

WAXSiS jobs generate graphical, data file, and structure file output. The results are summarized on a web page with a randomized URL to ensure user privacy. The web page shows the SWAXS curve, together with the fitted experimental curve (if present), as well as a Guinier fit. In addition, molecular representations of solute, hydration layer and excluded solvent are shown. The web page also presents a link to a gzipped tar ball that contains the following detailed results:

• Data and PDF files of the calculated SWAXS curve and Guinier fit, as well as the fitted experimental curve (if provided). Figure 3A shows a few examples of SWAXS curves computed by WAXSiS, together with fitted experimental curves. Experimental data shown in Figure 3A was taken from the Small Angle Scattering Biological Data Bank (26) (SASDA82 and SASDAK6) and from (12).

• A log file reporting details on the job options, ligand handling, MD simulation, SWAXS calculation, buffer subtraction scheme, fitting results and Guinier fit.

• A PDB file complete.pdb of the initial frame of the MD simulation system, containing the solute (including hydrogen atoms), water and counter ions. In addition, PDB files of the solute and hydration layer as well as of the excluded solvent are included, taken from the first frame of the MD trajectory after equilibration.

• A dummy topology file in Gromacs format. This file lists which element and hence, atomic form factor, was assigned to each atom. The atom numbers correspond to the atoms in complete.pdb. The user may consult this file to check if WAXSiS assigned the correct elements to uncommon atoms, such as exotic metal ions.

SWAXS curves from uploaded MD trajectories

As an additional feature, WAXSiS computes SWAXS curves based on explicit-solvent MD simulations conducted by the user. The user uploads the complete simulation system as a PDB file, the simulation trajectory (in Gromacs format), as well as an index file (in Gromacs format) defining the solute and the solvent. Details are provided in the Help section of WAXSiS. It is important to note that the simulation system should contain a rather large amount of water to ensure that the constructed envelope fits into the simulation cell.