Programming new geometry restraints: parallelity of atomic groups.

Improvements in structural biology methods, in particular crystallography and cryo-electron microscopy, have created an increased demand for the refinement of atomic models against low-resolution experimental data. One way to compensate for the lack of high-resolution experimental data is to use a priori information about model geometry that can be utilized in refinement in the form of stereochemical restraints or constraints. Here, the definition and calculation of the restraints that can be imposed on planar atomic groups, in particular the angle between such groups, are described. Detailed derivations of the restraint targets and their gradients are provided so that they can be readily implemented in other contexts. Practical implementations of the restraints, and of associated data structures, in the Computational Crystallography Toolbox (cctbx) are presented.

Introduction

The restrained refinement of atomic models is performed through minimization of a target function. First introduced in the macromolecular crystallography field by Konnert (1976 ▸) and Jack & Levitt (1978 ▸), this target is a weighted sum of several terms including describing the model fit to the experimental data:In the case of refinement of atomic coordinates, the terms describe the model geometry and are known as geometric or stereochemical restraints.

At low resolution there may be many models that may fit the experimental data equally well, but not all of these models may be chemically reasonable. This is why additional information needs to be used to restrict the models to those that conform to expectations about the correct chemical structure. These expectations may originate from various sources, such as small-molecule crystallographic studies, high-quality macromolecular structures or quantum chemistry calculations, forming the basis for libraries of a priori knowledge. At lower resolution, there are fewer experimental data available and therefore a larger number of restraints are required to ensure a high-quality refined model. Therefore, low-resolution structural studies, which have become more and more frequent, require additional restraints (see, for example, Smart et al., 2012 ▸; Headd et al., 2012 ▸; Brown et al., 2015 ▸) to complement the traditionally used restraints, such as bond lengths, bond angles, dihedral angles, non-bonded interactions, and planarity and chirality restraints.

As an example, nucleic acids have a number of base pairs, some of which are parallel or nearly parallel to each other. This knowledge may serve as a restraint that can be applied to the refinement of such models. Adding a new restraint into the refinement target means developing a corresponding algorithm to calculate its value, as well as an algorithm to calculate the gradient of this restraint function. We use the example of the parallelity restraint, implemented in cctbx (Grosse-Kunstleve et al., 2002 ▸), to illustrate the general approach of adding new restraints in cctbx and PHENIX (Adams et al., 2010 ▸) and also to reinforce the general principles of crystallographic refinement programs. This article describes the algorithms in a ready-to-program manner, with sufficient mathematical and computational details to serve a didactic purpose as well. Practical examples of the application of this new restraint function to model refinement will be presented elsewhere.

Each individual target in equation (1) is the sum of a large number of identical terms, each depending on a small number of parameters. For example, a diffraction target usually depends on a structure factor amplitude or the parameters of its phase distribution, a density target depends on the density values in a few grid points around a particular atom (see, for example, Afonine & Urzhumtsev, 2004 ▸, and references therein), a geometric restraint depends on the coordinates of a few atoms involved in particular restraint, etc. Such a composition of the target function (1) is very important for several reasons. First, it allows the manipulation of each kind of restraint, e.g. selection of the atoms involved or an update of the function. Second, each restraint is relatively simple and easy to program. Third, each restraint can be weighted individually. Finally, the gradient of the crystallographic target used for the minimization can be calculated as a sum of gradients from the individual terms.

Several general points are important to note when designing refinement programs (Lunin & Urzhumtsev, 1985 ▸):

(i) Refinement targets may be functions of the atomic parameters and values of electron density (or structure factors), in combination with the restraints and/or constraints imposed on the atomic model. Each component of the overall target is a relatively simple function when it is expressed through the parameters of the appropriate type.

(ii) These parameters are calculated from each other as a chain of consecutive transitions.

(iii) The gradient calculations are performed step by step using the chain rule, inverting the chain of steps used to calculate the target (Baur & Strassen, 1983 ▸; Kim et al., 1984 ▸). This procedure results in exact gradient values that require approximately the same amount of time to calculate as a single value of the target, and this is independent of the number of parameters to be refined.

(iv) The gradient calculation does not require an analytical function; it may be any algorithm that provides a value for the corresponding target. At each refinement step, the target or gradient calculation is a local operation. It takes the numerical values obtained at the previous step, calculates the result of the step and passes the numerical values to the next step. Such a scheme allows one to avoid complicated mathematical equations connecting the final target or gradient value with the initial model parameters.

In order to introduce a parallelity restraint into refinement, the considerations stated above require that we propose an algorithm to calculate a numerical value describing by how much two atomic groups that form pseudo-planes (approximately belong to two planes) are not parallel to each other. This algorithm should be numerically robust enough to handle a broad range of model parameters, from good (towards the end of refinement) to distorted (at the beginning of refinement).

Possible parallelity restraints and their gradient

Definition of the ‘best plane’

Let us consider two pseudo-planar atomic groups (groups that should be planar but may not currently be so), m = 1, 2, containing atoms defined by their Cartesian coordinates, r = (x , y , z ), k = 1, …, K, in the same coordinate system. For each such group we need to define the corresponding ‘best plane’. A plane is defined by a vector N normal to it and by a point C that belongs to it. For a given atomic group, the plane can be characterized by the sum of the squared deviations of the atoms from it (Urzhumtsev, 1991 ▸; Blanc & Paciorek, 2001 ▸; Brown et al., 2015 ▸), and these deviations may be weighted if desired. The ‘best plane’ is the one that corresponds to the minimum of this value (N and C are variables),where w is a weight (e.g. a unit value or an atomic mass) for the contribution of an atom with coordinates r to the target. Here and below, a centred dot ‘·’ represents the scalar (dot) product of two vectors, and is the normalized vector N .

Two facts are useful (see, for example, Urzhumtsev, 1991 ▸). First, the minimum of equation (2) is reached when C is the weighted geometric center of the groupSecond, if for a given atomic group we calculate the minimal eigenvalue of its inertia matrix (see §3.2), the corresponding eigenvector is normal to the best plane defined by condition (2) and therefore it can be taken as N .

The angle between two atomic groups can be described by the angle θ between two vectors normal to the best planes for the corresponding groups. This angle varies between the limits 0 ≤ θ ≤ π/2. This angle is expected to be zero if we require the planes to be parallel, which is the most frequent restraint of this kind. We therefore need to design a function which reaches its minimum when θ = 0. Note that there are multiple functions which could be constructed such that θ0 = 0 is a minimum. Also, the particular restraint function needs to describe the stereochemical knowledge near the point of interest (minimum) and have a desired behavior far from it. Additionally, the computation of the restraint and its derivatives must have no numerical irregularities. Below, we analyze some possible restraints that describe the angle between two planar groups, noting that different restraints can be used in different circumstances.

Analysis of various possible parallelity targets

The computationally simplest target to restrain the angle between normal vectors iswhere θ = θ(n 1, n 2) is the angle between two unit vectors n 1 and n 2, θ0 is the target angle, and w is the individual weight factor for a particular pair of planes. By default, phenix.refine (Afonine et al., 2012 ▸) uses the form w = 1/σ2, where σ is the root mean-square deviation from the ideal value obtained by some means. Introducing the angle θ0 in equation (4) will force the planes to form the prescribed angle between them. For example, as suggested by Brown et al. (2015 ▸), a restraint with θ0 = π/2 is important when working with aromatic groups. For simplicity, we continue to call this restraint ‘parallelity’, even though this is not completely appropriate. We note that cosθ = n 1·n 2 and sinθ = ||n 1 × n 2||.

For small deviations from θ0, restraint (4) behaves as the quadratic function w(θ − θ0)2. However, using the quadratic function of the angle itself [see in particular equation (14) of Brown et al. (2015 ▸)] requires its derivatives, i.e. the derivative of arccos(n 1·n 2), which is undefined for θ = θ0. Also, its calculation is computationally unstable for θ ≃ θ0. For this reason, we have excluded this kind of target from our consideration.

Also, a seemingly plausible restraint function that one may wish to avoid isThis function is easy to calculate but it has very different behavior across different ranges of the parameter θ0. For example, it is similar to a second-order polynomial near θ0 = π/2 and to a fourth-order polynomial near θ0 = 0. Such heterogeneity is generally disadvantageous for minimization.

Another known target for the parallelity restraint is[equation (13) of Brown et al. (2015 ▸); the notation of the present article is used]. It imposes a rigid condition on the parallelity of two planes, n 2 = n 1, simultaneously with planarity of the corresponding groups. This constraint-type condition may be useful, for example, to simplify the parameterization of the model, but in practice nucleic acid bases may bend from one base pair to another, or they may be twisted inside a given base pair. As a consequence, it is more convenient to have planarity and parallelity restraints defined separately, allowing for more flexibility if desired, in particular allowing the interplanar angle to vary from one planar group to another.

When the function is far from the minimum (e.g. when the current structure is distorted with respect to the restraints), one may wish that the target be relatively insensitive to the variation in the angle θ. This can be achieved by using various kinds of ‘top-out’ functions (Dennis & Welsch, 1987 ▸; Murshudov et al., 2011 ▸; Smart et al., 2012 ▸; Headd et al., 2012 ▸). A direct application of this idea results inwhere Ω is a constant (Fig. 1 ▸). Using (θ − θ0)2 or (cosθ − cosθ0)2 in equation (7) instead of cos(θ − θ0) − 1 is not desirable for the reasons discussed above. In fact, we can use directly a sigmoid form of the cosine function and introduce the ‘top-out type’ targets in a simpler form:orwhere n > 2. These functions (Fig. 1 ▸) have harmonic behavior when θ ≃ θ0, while for large values of the argument they are constant [equation (9)] or nearly constant [equation (8)].

Figure 1

Examples of various parallelity targets, θ0 = 0. Heavy line with no markers: f(θ) = 1 − cosθ, equation (4). Line with filled triangles: f(θ) = 1 − cos2θ, equation (8). Line with open triangles: f(θ) = 1 − cos4θ for θ ≤ π/4, equation (9). Lines with black (Ω = 1) and open (Ω = 2) circles: f(θ) = Ω2{1 − exp[(cosθ − 1)/Ω2]}, equation (7). Line with filled squares: f(θ) = (1 − cosθ)4, equation (10). Line with asterisks: f(θ), slack function, equation (11).

A further useful modification of target function (4) can be achieved by introducing the ‘slack’ parameter (square well-like potential), which scores small deviations from the target value equally. This can be achieved in a ‘soft’ way (Fig. 1 ▸) as or following explicitly the slack definition (Fig. 1 ▸) as

Practical calculation of the parallelity targets and their gradient

All the restraint functions mentioned above are similar in the sense that they require calculation of the normal vectors to the optimal planes of the atomic groups; then their dot product and finally a vector product of these vectors are calculated to obtain the target value.

A number of other restraints, such as group planarity or the distance between parallel groups, are naturally expressed through the same types of parameters used for the parallelity restraint (Appendix A ). These parameters are the center of the atomic group, and the eigenvalues and eigenvectors of its inertia matrix (generally speaking, all three eigenvalues and all three eigenvectors may be required); a particular target may use only some of the parameters mentioned. Thus, the calculation of these parameters from the initial set of atomic coordinates is common to all such restraints on planar groups (Fig. 2 ▸). For this reason, in §3 we provide details of these common steps.

Figure 2

The overall calculation scheme. (a) The steps used to calculate parameters common to all targets. The letters a–f from top to bottom refer to the steps for calculating the target described in §3. The same steps in the direction from bottom to top refer to the steps for calculating the derivatives described in §4. The steps are applied independently to each atomic group required to be planar. (b) The particular parallelity targets that can be constructed on the basis of the calculated parameters of an atomic group.

Assuming we have a function that we can calculate numerically, the principle of efficient gradient calculation (Baur & Strassen, 1983 ▸; Kim et al., 1984 ▸) suggests inverting the steps of the function’s calculation and obtaining the gradient step by step using the chain rule. This would result in the exact value of the gradient with respect to the original parameters, and this calculation will take no more than four times the CPU time required to calculate a single value of the function. This factor of four is invariant with respect to the number of parameters and the function type. In crystallographic practice, this factor is typically closer to one than to four (Lunin & Urzhumtsev, 1985 ▸). §4 illustrates this principle by calculating the gradient of the parallelity restraints. It is easy to see that the calculation of the gradient of other restraints for pseudo-plane groups differs only in the first step.

Planarity targets calculation

This section describes the calculation of the target value for the suggested restraint. Fig. 2 ▸ shows the overall scheme, where the letters a to f refer to the corresponding steps described below.

Centered coordinate system (Step a)

We start by calculating the atomic positions in the coordinate system with the origin shifted into the geometric center C of the group, equation (3). This operation is similar for all pseudo-plane groups, and in what follows we omit the group index m to simplify the formulae. The new parameters are the Cartesian coordinates (X , Y, Z) of the atomic positions: and the Cartesian coordinates (C, C, C) of the weighted geometric center of the group.

Inertia matrix (Step b)

For a given atomic group, its inertia matrix isThis matrix is symmetric, S = S, S = S, S = S, but for convenience in subsequent calculations we refer to all nine of its elements. Matrix (13) is positive semidefinite, having always three non-negative eigenvalues, 0 ≤ λ+ ≤ λ− ≤ λ0, with the corresponding eigenvectors orthogonal to each other. The result of this step contains nine elements of matrix (13) and three coordinates of the center [equation (12)] taken from the previous step.

It is important to remember that at each step we calculate the values of the new parameters and do not keep their analytical expression.

Characteristic equation (Step c)

The eigenvalues 0 ≤ λ+ ≤ λ− ≤ λ0 of the inertia matrix (13) are the roots of the cubic characteristic equationwith the coefficients Further steps of the procedure use both the coefficients of the cubic equation (14) and the elements of the inertia matrix (13). To avoid confusion when calculating the gradients [see, for example, equation (18) below], it is more convenient to introduce new variables formally:Thus, the new set of parameters contains 15 values: three coordinates (C, C, C) of the center of the atomic group that are simply taken from the previous step, three equation coefficients a, b, c, and nine elements of the inertia matrix .

Eigenvalues of the inertia matrix (Step d)

Different algorithms may be applied to calculate the eigenvalues and choose the minimum of them; we use the one described by Urzhumtsev (1991 ▸), as it gives analytically the minimal real eigenvalue that is required for the planarity and parallelity restraints. Appendix B gives an extended description in the notation of the present paper. The coordinates of the center of the atomic group [equation (12)] and the matrix elements [equation (16)] are not used at this step and pass to the next steps without change.

Eigenvectors of the inertia matrix (Step e)

Different algorithms may be applied to calculate the eigenvectors N of the inertia matrix that are solutions of the matrix equationwith the eigenvalues λ obtained in the previous step. Appendix B gives technical details of the method used in our program. For an atomic group that is approximately planar, λ+ < λ−, the solution of equation (17) with λ = λ+ defines unambiguously a single direction N + of the normal to the best plane.

Note that other restraints on the plane groups may require knowledge of other eigenvalues and/or eigenvectors (Appendix A ). Therefore, in the general case, the result of this step is a set of parameters fully describing a pseudo-planar atomic group, namely the three coordinates (C, C, C) of its center, all three eigenvalues 0 ≤ λ+ ≤ λ− ≤ λ0 and all three eigenvectors N +, N −, N 0 of the inertia matrix.

Parallelity targets (Step f)

To calculate a planarity restraint for two pseudo-planar atomic groups, each described by the Cartesian coordinates of the vector N 1 = N 1+ and N 2 = N 2+ normal to their respective best planes, we may need to follow several intermediate steps (Appendix D ). We start (Step f1) from normalization of the normal vectors, {N 1, N 2} {n 1, n 2}. This normalization includes choosing the correct direction of the normalized vectors so that the value of the angle θ between the normal vectors is not larger than π/2. Once the normalized vectors are known, we calculate (Step f2) the cosine and sine of the angle θ between them, which we denote C θ and S θ, respectively. Finally (Step f3), we calculate the value of the chosen target.

Gradient calculation

Gradient with respect to the plane parameters

We start by inverting the formulae to calculate the chosen target (Step f3 in Fig. 2 ▸) and obtain and (Appendix D ).

Inverting the previous step (Step f2), we take as input the values of and (and not their analytical expressions), whatever the target R parallelity is. As a result, we obtain (, , ) and (, , ). Finally, inverting the normalization (Step f1) we obtain (, , ) and (, , ).

Working with other restraints on plane groups, we can apply the same kind of procedure as Steps f3 to f1 above. For all these targets, we arrive at a set of partial derivatives of a target with respect to the following:

(i) The three eigenvalues for each pseudo-planar group.

(ii) The coordinates of the three eigenvectors for each pseudo-planar group.

(iii) The coordinates of the center of each pseudo-planar group.

For the parallelity restraints discussed above, all these derivatives are equal to zero except those for the coordinates of N 1+ and N 2+. Further steps of the gradient calculation use these partial derivatives as input data whatever the target is. These steps are common to all restraints and all atomic groups. Therefore, starting with the next section we omit the index used to indicate an atomic group and the index of a particular plane-group restraint for which all the partial derivatives above have been calculated. Partial derivatives with respect to the coordinates (C, C, C) will not change until the last step.

Gradient with respect to the coefficients of the characteristic equation (Steps e and d)

As mentioned in §§3.4 and 3.5, different algorithms exist for calculating the eigenvalues and eigenvectors of the inertia matrix from the coefficients of the characteristic equation (14). According to the main scheme, we invert these algorithms to calculate the corresponding derivatives. Appendix C provides these for the particular method of calculation chosen in our program. The output of this step contains the partial derivatives of the target R with respect to the three coefficients a, b, c of equation (14), to the nine elements , …, of equation (16), and to the three coordinates (C, C, C) of the center of the atomic group that are taken, with no change, from the previous step.

Gradient with respect to the elements of the inertia matrix (Step c)

Using the partial derivatives of a, b, c [equation (50) in Appendix C ] and of the elements [equation (43) in Appendix C ] with respect to the elements of the matrix S [equation (13)], we obtainand similar expressions for the other six derivatives easily obtained from equation (18) by the corresponding cyclic substitution of the indices, e.g. to obtain we make the following substitutions in the second equation: , , , , .

Gradient with respect to the coordinates in the centered system (Step b)

We obtain the derivatives with respect to the atomic coordinates, k = 1, …, K,

Gradient of restraint with respect to the original coordinates (Step a)

Using the partial derivatives of equations (3) and (12) and the notation W = ,and knowing that all cross-term derivatives are equal to 0, we obtain the required gradient with respect to the original atomic coordinates

Implementation of parallelity restraints in cctbx

Computational Crystallography Toolbox

The Computational Crystallography Toolbox (cctbx) is an open-source library with numerous functionalities necessary to construct crystallography-oriented programs. In fact, it finds applications beyond crystallography, such as small-angle scattering (sastbx; Liu et al., 2012 ▸) or cryo-EM (Afonine et al., 2013 ▸). cctbx contains tools to handle various objects such as crystal symmetry, diffraction data, maps (crystallographic or cryo-EM) and atomic models, and provides an array of building blocks necessary to construct complex procedures such as model refinement (Afonine et al., 2012 ▸) or diffraction data analysis (Sauter et al., 2013 ▸; Parkhurst et al., 2014 ▸). In particular, to support macromolecular refinement algorithms cctbx includes tools to handle various geometry restraints (Grosse-Kunstleve & Adams, 2004 ▸). The main set of geometry restraints includes those on covalent bonds and angles, chirality, dihedral (torsion) angles, non-bonded interactions, planar groups, and now parallelity. The overall restraint target, the sum in equation (1), is the sum of contributions from all restraint typesIf additional information about the molecule or its environment is available and deemed necessary to use, additional terms can be added to equation (22) (Echols et al., 2010 ▸; Headd et al., 2012 ▸, 2014 ▸). Such additional information may be the regular conformations of a protein’s secondary structure, non-crystallographic symmetry or similarity to a previously solved structure. Technically, adding additional restraints means adding additional terms to equation (22), which in turn means that an algorithm for calculating the additional restraint target and its gradient with respect to atomic coordinates, and the data structures for holding information about the atoms involved, are required.

The implementation of the parallelity restraint in cctbx is very similar to other restraints (angle, dihedral angle, planarity etc.) and therefore may serve as an example of the common framework for restraints.

General approach to defining stereochemical restraints in cctbx

The source code for commonly used restraints is located at http://sourceforge.net/p/cctbx/code/HEAD/tree/trunk/cctbx/geometry_restraints in the files bond.h, bond_bpl.cpp, angle.h, angle_bpl.cpp, dihedral.h, dihedral_bpl.cpp, chirality.h, chirality_bpl.cpp, nonbonded.h, nonbonded_bpl.cpp, planarity.h, planarity_bpl.cpp, parallelity.h and parallelity_bpl.cpp.

Files with the name _bpl.cpp contain boost.python (Abrahams & Grosse-Kunstleve, 2003 ▸) descriptions of classes and functions for the restraint, and the bpl in the file names stands for ‘boost.python library’. Generally, these files contain boost.python wrappers for classes describing the restraint and the associated data structure (termed proxy; see below) holding the information necessary to calculate the restraint. Files with the name .h contain the actual implementation of classes and functions for the particular type of restraint.

Proxy class

The calculation of a particular restraint requires information about the atoms involved, such as the atomic coordinates, and the restraint target and weight. Since the restraints are built once at the beginning of refinement, while the atomic coordinates change during refinement, information about the atoms that participate in a restraint and the actual atomic coordinates is decoupled. The data structure that holds the restraint information is called a proxy, and it contains integer indices that can be used to access a particular atom in an array of atoms and therefore individual coordinates, along with some other restraint-specific information.

For example, the bond proxy contains a pair of integer indices (atom numbers in an atom list), the restraint target value and weight, the slack and limit values, a flag that specifies which function to use (least-squares or top-out) and a rotation matrix operator in case bonded atoms belong to different crystallographic symmetry copies, and an indicator to distinguish the origin of a particular restraint (e.g. covalent or hydrogen bond). Since each proxy object contains all the parameters necessary to define a particular restraint, each restraint can be parameterized individually.

The parallelity proxy is organized in a very similar manner to the planarity proxy, the difference being that the atom numbers are split and contained in two groups because they define two separate planes with an arbitrary number of atoms in each one.

Usually, one constructs multiple restraints of one type for a structure, so it is convenient to have a structure to store an array of proxies of one type and perform various operations on them. Similarly to other restraint types, for the parallelity restraints there is an array that can contain multiple proxy objects, called the shared_parallelity_proxy. Operations on proxy arrays include selection (this allows selection of a subset of restraints corresponding to a part of the structure), deletion (restraints can be deleted for a part of structure), and calulation of root mean-square deviations of the model from target values, absolute differences between the restraint target and actual model values, the restraint target value, gradients and so on.

Class to handle parallelity restraint

The parallelity restraint itself [equation (4)] and its modifications [equations (7) and (11)] are implemented as a C++ class parallelity available in Python via boost.python from cctbx.geometry_restraints.parallelity. It has two constructors. One requires providing all parameters to define a parallelity proxy (such as target value, type of function etc.) and two sets of Cartesian atomic coordinates corresponding to the two planes. The other constructor takes an array with Cartesian coordinates of atoms and an instance of the parallelity proxy class that contains the rest of the required information.

All constructors contain a call to the init_deltas() function which calculates all necessary values, such as the overall restraint value and the gradients. This means that instantiation of a parallelity object leads to the immediate calculation of the target value and gradients for the restraint.

To simplify the handling of the slack functionality and to make it universal for harmonic and top-out potentials, the conditions of equation (11) are evaluated first and θ0 adjusted accordingly before any other calculations. This makes it possible to avoid the gradient calculation if it is not necessary, for example when |θ − θ0| ≤ slack, and also to use all of the derivative equations for harmonic and top-out potentials without change.

In addition to the two constructors, the parallelity class has member functions residual() and gradients() that return the value of the restraint function and gradients with respect to the atomic coordinates. Since all necessary calculations are already performed in the constructor, these functions return precalculated values.

Discussion

Applying the general principles of refinement program construction, it is straightforward to build various stereochemical restraints. However, a careful analysis is required to choose an appropriate function for the restraint target. This choice is based on a function’s ability to provide the desired behavior, both near the minimum and far from it. Possible singularities when calculating the target and its derivatives must also be taken into account. With a view to improving the refinement of atomic models, especially at low resolution, we have developed new restraints to control the angle between two approximately planar atomic groups. As a particular and more frequent case, this restraint imposes the groups to be parallel or near parallel. We have proposed a practical algorithm to calculate this target and its gradient. This algorithm has been implemented in the cctbx open-source library, making it available in PHENIX and to the broader crystallographic and macromolecular structure community. All calculation details are provided, making it easy to add other types of restraint in the future.